text:
How to convert fasta file to tab delimited file
I have a fasta file like
>sample 1 gene 1
atgc
>sample 1 gene 2
atgc
>sample 2 gene 1
atgc

I want to get the following output, with one break between the header and the sequence.
>sample 1 gene 1   atgc
>sample 1 gene 2   atgc
>sample 2 gene 1   atgc

Thanks so much everyone. You're making it hard to choose. I wanted something for a multi line fasta so both terdon and Chris scripts are correct. So I will defer to the saying first come first served.

If you have multi-line fasta files, as is very common, you can use these scripts1 to convert between fasta and tbl (sequence_name <TAB> sequence) format:

FastaToTbl
#!/usr/bin/awk -f
{
        if (substr($1,1,1)==">")
      if (NR>1)
                     printf "\n%s\t", substr($0,2,length($0)-1)
      else
       printf "%s\t", substr($0,2,length($0)-1)
            else
                    printf "%s", $0
}END{printf "\n"}

TblToFasta
#! /usr/bin/awk -f
{
  sequence=$NF

  ls = length(sequence)
  is = 1
  fld  = 1
  while (fld < NF)
  {
     if (fld == 1){printf ">"}
     printf "%s " , $fld

     if (fld == NF-1)
      {
        printf "\n"
      }
      fld = fld+1
  }
  while (is <= ls)
  {
    printf "%s\n", substr(sequence,is,60)
    is=is+60
  }
}

Save those in your $PATH, make them executable, and you can then do:
$ cat file.fa
>sequence1
ATGCGGAGCTTAGATTCTCGAGATCTCGATATCGCGCTTATAAAAGGCCCGGATTAGGGC
TAGCTAGATATCGCGATAGCTAGGGATATCGAGATGCGATACG
>sequence2
GTACTCGATACGCTACGCGATATTGCGCGATACGCATAGCTAACGATCGACTAGTGATGC
ATAGAGCTAGATCAGCTACGATAGCATCGATCGACTACGATCAGCATCAC
$ FastaToTbl file.fa
sequence1   ATGCGGAGCTTAGATTCTCGAGATCTCGATATCGCGCTTATAAAAGGCCCGGATTAGGGCTAGCTAGATATCGCGATAGCTAGGGATATCGAGATGCGATACG
sequence2   GTACTCGATACGCTACGCGATATTGCGCGATACGCATAGCTAACGATCGACTAGTGATGCATAGAGCTAGATCAGCTACGATAGCATCGATCGACTACGATCAGCATCAC

And, to get the Fasta back:
$ FastaToTbl file.fa | TblToFasta
>sequence1
ATGCGGAGCTTAGATTCTCGAGATCTCGATATCGCGCTTATAAAAGGCCCGGATTAGGGC
TAGCTAGATATCGCGATAGCTAGGGATATCGAGATGCGATACG
>sequence2
GTACTCGATACGCTACGCGATATTGCGCGATACGCATAGCTAACGATCGACTAGTGATGC
ATAGAGCTAGATCAGCTACGATAGCATCGATCGACTACGATCAGCATCAC

This can be a very useful trick when searching a fasta file for a string:
TblToFasta file.fa | grep 'foo' | FastaToTbl

If you really want to keep the leading > of the header (which doesn't seem very useful), you could do something like this:
$ perl -0pe 's/\n//g; s/.>/\n>/g; s/$/\n/;' file.fa
>sequence1 ATGCGGAGCTTAGATTCTCGAGATCTCGATATCGCGCTTATAAAAGGCCCGGATTAGGGCTAGCTAGATATCGCGATAGCTAGGGATATCGAGATGCGATAC
>sequence2 GTACTCGATACGCTACGCGATATTGCGCGATACGCATAGCTAACGATCGACTAGTGATGCATAGAGCTAGATCAGCTACGATAGCATCGATCGACTACGATCAGCATCAC

But that will read the entire file into memory. If that's an issue, add an empty line between each fasta record, and then use perl's paragraph mode to process each "paragraph" (sequence) at a time:
perl -pe  's/>/\n>/' file.fa | perl -00pe 's/\n//g; s/.>/\n>/g; s/$/\n/;'

1Credit to Josep Abril who wrote these scripts more than a decade ago.

I added an answer just to show possible improvements to the 2 awk scripts which IMHO are over-complicated with too much code duplication and other issues.

There is a very simple BioPython solution, that is minimal, readable, and handles multi-line fasta:
from Bio import SeqIO

for record in SeqIO.parse('example.fa', 'fasta'):
    print('>{}\t{}'.format(record.description, record.seq))

assuming there is only one sequence line per record, use paste with two 'stdin'
cat your.fasta | paste - -

Note that this will fail if you have multi line sequences (as Pierre pointed out), but also if you have any blank lines in the file. You might also want to remove the UuOC: paste - - < file.fa.

You can use these commands:
perl -pe 's/>(.*)/>\1\t/g; s/\n//g; s/>/\n>/g' sequences.fa | grep -v '^$'

Explanation:

Append a tab to every header line
Join all lines
Split the single obtained line by the '>' character
Remove the empty line (the first line is empty due to the fact that '>' is the first character of the FASTA file)

A very useful tool for this kind of data manipulation is bioawk:
$ bioawk -c fastx '{print ">"$name" "$comment"\t"$seq}' test.fa
>sample 1 gene 1    atgc
>sample 1 gene 2    atgc
>sample 2 gene 1    atgc

bioawk is based on awk, with added parsing capabilities. Here, we tell that the format is fasta or fastq with -c fastx, and this makes the $name (between ">" and the first blank character), $comment (after the first blank character) and $seq (the sequence, in one line) variables available within awk instructions.
See for instance this answer for another use case.

Where possible, I recommend using a dedicated parsing library, rather than hacking a parser together: as you can see in the other answers, parsing even simple formats gets complex pretty quickly if you value correctness.
Here’s a small R script that does what we need, using ‘seqinr’:
#!/usr/bin/env Rscript
suppressPackageStartupMessages(library(seqinr))
parsed = read.fasta(file('stdin'), as.string = TRUE)
table = data.frame(unlist(parsed), row.names = sapply(parsed, attr, 'Annot'))
write.table(table, stdout(), sep = '\t', quote = FALSE, col.names = FALSE)

Save it as fasta-to-tsv, make it executable, and use it as follows:
fasta-to-tsv < input.fasta > output.tsv

Equivalent code of similar length can be written in Python or Perl.

Could you also explain what packages need to be installed in order for R to do this? seqinr isn't part of vanilla R.

@terdon I’m slightly confused: ‘seqinr’ is the package name, so that’s the one to install.

Well, sometimes the name of the package isn't the same as the library you load from it. And in any case, R being the horribly complicated mess that it is, you can never know (or I can't, anyway) whether the package can be installed with install.package or needs to be done through bioconductor or whatever else may be the case. R package management is very far from trivial, so I thought it would be useful to explain how to install the packages necessary to run your code. I've edited the command in now. I'd upvote, but I already did so back in 2017 :)

@terdon “sometimes the name of the package isn't the same as the library you load from it” — No, for R it always is, by definition. As for installation, I’ll give you that there’s a special case for packages not on CRAN but it only really makes sense to spell out the exceptions, not the rule. You wouldn’t mention installation instructions for Perl packages on CPAN, nor for Python on PyPI, NPM for JavaScript or Rust crates. I’ve undone your edit because (a) it’s redundant, and (b) it was wrong. Never install packages via sudo R unless you’re the sysadmin maintaining a cluster.

Ah, I tried it without but it failed when trying to install to my user's local library, so I assumed the package could only be installed system wide. But that just proves why it is helpful to include the installation command. I don't know R (and hate it with a passion) so all these things that you mention which are obvious to you, are not to me. And possibly others, which is why yes, I always include commands to install packages in my answers.

This could be easily done by seqkit fx2tab
seqkit fx2tab seq.fa

However, seqkit will not print the "greater than" symbol (">"). If you do need the symbol:
seqkit fx2tab seq.fa | sed 's/^/>/g'

Remove empty records (description without sequence):
awk '$2{print RS}$2' FS='\n' RS=\> ORS= f1.fa > f2.fa

Remove blank lines:
sed '/^$/d' f2.fa > f3.fa

Convert multi-line fasta to single-line fasta:
awk '/^>/ {printf("\n%s\n",$0);next; } { printf("%s",$0);} END {printf("\n");}' f3.fa > f4.fa

Finally, @Pierre solution:
cat f4.fa | paste - - > f.txt

In cases where there is no sequence wrapping and each sequence occupies only a single line, the following shell command is probably going to be fastest, easiest, and most convenient.
paste - - < your.fasta > your.new.fasta

Oops, looks like Pierre already suggested this. :-)

You don't need to leave duplicate answers lying around, you could delete it.

Just suggested improvements to the awk scripts in @terdon's answer, using any POSIX awk:
$ cat FastaToTbl
#!/usr/bin/env bash

awk -v OFS='\t' '
    {
        if ( /^>/ ) { out = (NR>1 ? ORS : "") substr($0,2) OFS }
        else        { out = $0 }
        printf "%s", out
    }
    END { print "" }
' "${@:--}"

$ cat TblToFasta
#!/usr/bin/env bash

awk -F'\t' '
    {
        gsub(/.{60}/,"&"ORS,$2)
        sub(ORS"$","",$2)
        print ">" $1 ORS $2
    }
' "${@:--}"

$ ./FastaToTbl file.fa
sequence1       ATGCGGAGCTTAGATTCTCGAGATCTCGATATCGCGCTTATAAAAGGCCCGGATTAGGGCTAGCTAGATATCGCGATAGCTAGGGATATCGAGATGCGATACG
sequence2       GTACTCGATACGCTACGCGATATTGCGCGATACGCATAGCTAACGATCGACTAGTGATGCATAGAGCTAGATCAGCTACGATAGCATCGATCGACTACGATCAGCATCAC

$ ./FastaToTbl file.fa | ./TblToFasta
>sequence1
ATGCGGAGCTTAGATTCTCGAGATCTCGATATCGCGCTTATAAAAGGCCCGGATTAGGGC
TAGCTAGATATCGCGATAGCTAGGGATATCGAGATGCGATACG
>sequence2
GTACTCGATACGCTACGCGATATTGCGCGATACGCATAGCTAACGATCGACTAGTGATGC
ATAGAGCTAGATCAGCTACGATAGCATCGATCGACTACGATCAGCATCAC

FastaToTbl will actually work in any awk but TblToFasta requires a POSIX awk for support of regexp intervals like {60}. If you have a very old, pre-POSIX awk that doesn't support regexp intervals then get a new awk but if that's impossible for some reason then change TblToFasta to the following and then it'll also work in any awk:
$ cat TblToFasta
#!/usr/bin/env bash

awk -F'\t' '
    BEGIN {
        dots = sprintf("%"60"s","")
        gsub(/ /,".",dots)
    }
    {
        gsub(dots,"&"ORS,$2)
        sub(ORS"$","",$2)
        print ">" $1 ORS $2
    }
' "${@:--}"

Heh, I have been using these two scripts for decades, long before I knew the first thing about awk, and I've never bothered to go and review them. Your versions are much simpler, thanks!

It's hard to figure out why that second one especially was so complicated. I feel like it maybe was intended to handle additional types of input, not just what's output by the first script, e.g. maybe lines which have multiple chunks of sequence string separated by spaces? That loop on field numbers but doing something specific for field numbers 1 and NF-1 (which could also be 1?) and printing every other field, then falling into a second loop that's again printing substrings of the last field is baffling.

Ohhhhhh hang on - I get it. That first loop is trying to handle a sequence name that can contain spaces because they didn't set FS to \t. There's obviously far more concise ways they could have done that (but if we just set FS to \t then we don't have to). I'm assuming sequence names can't contain tabs.

This is an old post I have noticed and there are many offered solutions. Since it’s a frequently asked question, I thought it’s worth for me to mention that there is an overlooked tool set which contains a stand-alone program called faToTab in addition to many other useful bioinformatics tools.
faToTab inputFile.fasta outFileFasta_tab.txt
It’s a gold-chest in my opinion. Here are the links to utilities folder and details:
Description and Download instructions -
Binaries by machine type -
Link to the github page.
Anaconda installation is:
conda install -c bioconda ucsc-fatotab
conda install -c bioconda/label/cf201901 ucsc-fatotab

Thank you. I did not know there was a conda installation. It is convenient for sure.

In python I would do:
#...Suppose you have header information and
#...sequences stored in lists

header_info1 = [elements]
header_info2 = [elements]
sequences = [sequences]
index = enumerate(sequences)
table = open(pathtoyourfile.tsv,'w+')
for h1,h2,s,i in zip(header_info1,header_info2,sequences,index):
    table.write(f">{h1}\t{h2}\t{s}\n")
    if i+1==len(sequences):
       table.write(f">{h1}\t{h2}\t{s}")
table.close()

So basically I use f strings and I iterate over these three vectors that are of the same length. At the end of the iteration, I remove new line(\n) since it won't be needed anymore since you have not to write anything further.

How is this any different / better than the other pythonic answer? This looks a lot more verbose but doing actually less than the original one (assuming data is loaded already). It could be I am just missing something, but it would be good to explain why your solution might be favorable in some cases.

It is not that my solution is favorable in comparison to the other. As a developer, I found myself very often involved in situations where you have to bend the codes to my needs. Rather than understanding pre-made functions and losing time, I prefer to go on my own. Moreover, I saw that, commonly, f strings are not very well known yet, and in python-format function is used. I thought only it was a good occasion to share something that I commonly use.


text:
Find Patterns in Cluster
I have a heatmap and I would like to find some rectangles.

I have already used clustermap. But here, I can not calculate these rectangles.
The order of the data should not be changed.
This Code is not working:
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
tips2 = pd.read_csv('heatmapExport.csv')
sns.set(color_codes=True)
g = sns.clustermap(tips2, cbar="false")
plt.show()

Does somebody have an idea?
Thanks a lot
Here is an example of what I want:

What is pd, pandas? And sns? What are the reasons to select those rectangles and not others (without the reasons it will be hard to find a solution)? Why do you say the code is not working, how does it fail (which error, what is the output compared to the expected output...) ?

Yes, pd is pandas. I want to calculate the sum of the rectangles. I have a 2d array with float values. In this array I would like to find clusters. Optical are these clusters to be recognized in the heatmap. It would be good now if I had an algorithm that recognizes these figures and calculates the sum in the clusters.

It is therefore important that the points are not viewed independently of each other.

Sorry, but I don't understand. Do you want an algorithm that for the manually selected rectangles sum the values? Why don't you calculate it manually selecting the positions of the rectangles? BTW, what is the purpose of this/the biological question ?

Okay, I think I've said this wrong. Step 1: The algorithm is to form clusters. Step 2: The recheck should be entered into the heatmap. Step 3: Calculate the sum of the rectangles.

What is "the recheck"?  How do you define the rectangles?A rectangle is a cluster of samples and genes?

Oh sorry, recheck = rectangle... Sorry typo!

A rectangle is to form the cluster in which n * m data points are located

Can you give us a Minimal, Complete, and Verifiable example with some sample data, the rules to find the cluster ('recognizing optically' is not exact enough), and the expected results?

See the edited question above

Your example says "Cluster (for example)". Did you have criteria in mind to delimit the rectangles this way? Has the delimitation of rectangles an influence on the delimitation of others (for instance, by forbidding overlaps)? I think your question amounts to being able to formalize an intuition you have about how to form those rectangles.

Also, edit your question to add an explanation of what you mean by "This Code is not working". That's an important point.

I agree with @bli, If you could explain how do you want to make those clusters something could be done, but without knowing how do you want them we can't help you :(

I do not exactly know how to Form this clusters. This is one part of my Question. The criteria would be to find a Minimum number of clusters that "fit" the  I mean if you look at the heatmap. You can See this clusters without having any criteria. This is what i Want i do. But the Algorithm can Not Look at the points Independent of each other.

Does anyone have an idea how to deal with the problem?

Based on your description I think you should have a look at a technique called 'biclustering'.
The example on this page defines the goal of this technique as 'Finding subgroups of rows and columns which are as similar as possible to each other and as different as possible to the rest.'
Since your examples are python-based, you could check out scikit-learn's implementations of biclustering.

Many thanks for your response. I have implemented the algorithm as described. He does (almost) what I want. Now I have only one problem: I would not re-sort the rows and columns. For example, the algorithm says: Cluster 1 [line 1,3,4; Columns 1, 2, 3] Cluster 2 [line 2; Column 4]. I would like to know the result, however, the other way around. Lines 1.2 Columns 1 -> cluster 1; Lines 3.4 columns 2,3 -> clusters 2, ... etc. At all, it is important that the different Clusters Switch between rows and cols. If one Cluster has finished, no Datapoint should be add to this cluster

I would second scitest-learn as the approach

source: Stack Exchange

license: Creative Commons - Attribution Share-Alike - https://creativecommons.org/licenses/by-sa/3.0/

sort: votes


text:
Obtain CD (cluster of differentiation) marker expression information
There are 394 CD marker genes. Various resources (BD handbook, Abcam guide) contain information about which cell types express which CD markers. What is the best way to obtain or download this data so that it can be used with programs such as R?
How can we download similar annotation data on expression for any gene or cell type?

Most of the work of biologist is to characterize the cells and populations, mostly because they reason why they express/use them is mostly unkown. I am not aware of a single database listing markers of a cell line, but it would be quite useful.

I can offer this as a way to scrape the data from the Abcam resource:
library(rvest)
library(dplyr)

antigens <- read_html("http://www.abcam.com/primary-antibodies/human-cd-antigen-guide") %>%
  html_node("table") %>%
  html_table(header = TRUE)

# what targets expressed in "Langerhans" ?
antigens %>%
  filter(grepl("Langerhans", `Cellular Expression`)) %>%
  select(`Target Name`)

   Target Name
1         CD1a
2         CD1b
3         CD1c
4         CD1e
5         CD14
6         CD39
7         CD50
8         CD69
9         CD74
10        CD83
11       CD101
12       CD206
13       CD207

The data could be further processed by, for example, splitting the Cellular Expression column on either and or , (comma) to try and make one row per cell type.

all_licenses: Creative Commons - Attribution Share-Alike - https://creativecommons.org/licenses/by-sa/3.0/


text:
Simulate fast5 MinION run
I'm planning on using Oxford Nanopore MinION sequencer to test a script that I'm working on. However, I would not like to run multiple MinION runs to validate the tool. I'd like to simulate various minION runs and get raw reads and all pore events in a fast5 file.  Do you know of any Oxford Nanopore (MinION) sequencing simulator that allows to get fast5 as output?
I looked into NanoSim but it seems like there is no option to get fast5 as an output.

Possible duplicate of Tools for simulating Oxford Nanopore reads

Hi Andresito. This question seems to be a duplicated of the linked question by gringer. However you could [edit] it to ask about the fast5 output if you would like to have an answer (and to keep open the question).

Using a simulator for example MinION output doesn't make sense when there are an abundance of freely-available FASTQ and FAST5 sequences for general purpose use (and especially for the purposes of testing out algorithms). Simulation of sequences can only ever express the components of DNA that are explicitly modelled, which is never going to fully represent reality.
As long as you're doing your own private research, you can use the FAST5 sequences from Clive Brown's own genome:
https://github.com/nanoporetech/ONT-HG1
Bear in mind the license attached to that data:
https://github.com/nanoporetech/ONT-HG1/blob/master/LICENSE.md

License to use for research use only (i.e. non-commercial) governed by the law of England and Wales.
You agree to this License by downloading the data on this site whether in whole or in part.
You cannot re-distribute modified copies of these data, or derivative works, without my express written permission. Where such permission is granted, you must cite the original work and its Copyrights and you must include this license with any work that includes any of my genome.

For a more open license, you can use the nanopore-wgs reads here (derived from human cell line NA12878):
https://github.com/nanopore-wgs-consortium/NA12878

We encourage the reuse of this data in your own analysis and publications which is released under the Creative Commons CC-BY license. Therefore we would be grateful if you would cite the reference below if you do.

If you don't want human data, there's mouse parasite data here:
http://www.ebi.ac.uk/ena/data/view/PRJEB20824
And E. coli data here:
https://www.ebi.ac.uk/ena/data/view/PRJEB22772
... and if you just want 8h of sensor noise as a FAST5 file, I've got that here in one of my DropBox folders (as a ~70MB file):
https://www.dropbox.com/s/20j1anlz0motrs4/MINICOL233_20170118_FN__MN15625_sequencing_throughput_ONLL3127_84599_ch320_read10302_strand.fast5?dl=0
Reminder: if you do ever present the results of analyses based on a data source, it's a good idea to acknowledge the source (even if not expected). If nothing else, it will let the creator know that someone cares about their data, and perhaps encourage them to release more in the future.

Thank you very much @gringer. I will revise the provided sources and acknowledge accordingly. You are great!

I wrote a tool that will rerun an experiment. So all you would need to do is just download an ONT dataset from ENA or something and you can just rerun it as many times as you want.
https://github.com/mbhall88/taeper


text:
How can I compare structures of proteins responsible for dehydration tolerance in plants?
For my undergraduate research project, I'm trying to compare the structures of some proteins, specifically transcription factors crucial for dehydration tolerance in plants, in the hope of finding differences that make each protein more or less dehydration tolerant than the others. i want to compare the selected transcription factors at their structural level. structure of one transcription factor from different plant species such as rice,sorghum,maize etc. are considered. ultimately i'm hoping to find a relationship between these structures of proteins and their respective level of dehydration tolerance in each plant.
Is there any database or a site where I can find information regarding this?

Is it the proteins that are dehydration tolerant or the plants? Are you looking for a database of protein or plant features? And are you interested in the transcription factors or in the protein products of the genes regulated by the transcription factors? Please [edit] your question and explain what you are trying to do in more detail and tell us what you have so far.

i'm not looking for a database of plant features. hoping to find more information regarding the comparison of protein structures. thank you!

So you want to compare the structures of the transcription factors? Do you have a list of proteins of interest? Have you checked if they have known structures? Will you need to predict the structures? Have you looked at any structural alignment tools? More importantly, what makes you think the structures of the transcription factors will be relevant here? This sounds very much like an XY problem.

This project is done both in the protein sequence level and the structural level. I'm curious whether the plants with different levels of dehydration tolerance is maybe caused by the differences of responsible proteins at their structural level

OK, but you haven't answered my questions and we need the answers to be able to help. Do you have a specific list of proteins of interest and do you know if those proteins have known structures? And, again, it doesn't make much sense to be looking at the structure of the TFs, you should instead look at the structures of the protein products of the genes controlled by the TFs. The structure of the TFs will represent their function as transcription factors.

@terdon I believe your last comment is not entirely correct, if the TF changes conformation upon drought it would be extremely important to take the structure into account.

Have you explored in PDB?. If you already have a starting sequence, you can use it as input to search with phmmer using PDB as database. Once you have a curated set to analyse, you can visualize and align the structures with VMD, for example.


text:
How to provide a GFF/GTF file to Chromonomer
I am trying to use Chromonomer. It should be possible (according to the software page and a paper) to add a GFF/GTF file to the analysis.
Alas, neither the built-in help nor some light source code browsing were enough to find a way of passing the file to the software.
How can I use a GFF or a GTF file with Chromonomer?

I think the part of the code that should allow you to do this is in chromomer.cc.
Lines 2857-2983 form the function parse_command_line. Nowhere is there any reference to anything resembling parsing a GFF file, so I think the documentation is not in line with the code at this point.
For reference, I looked at version 1.07. It would be helpful if the code were on github.

The distribution of Chromonomer comes with a script for translation.


text:
Is there a database of disease prevalence?
I want to collect data on how common different diseases are in different populations (or, at least, globally). Ideally, it would give me a way of querying the database with a disease name, and would return the number of cases per N population.
My main objective is to have a general idea of how prevalent a disease is, I am not trying to track current outbreaks, so I am not too bothered if the information is a little out of date, as long as the database is still maintained and not abandoned.
I have found FIDD, but it seems to be no longer maintained. Are there any actively maintained, similar tools that have the information I am after?

For rare diseases, you can use Orphanet. You can download an xml file with prevalence from orphadata.
There is an OWL version of Orphanet called ORDO, you can browse it on the EBI Ontology Lookup Service. See for example Arrhythmogenic right ventricular cardiomyopathy and point prevalence relationships in the lower right box. There are a variety of programmatic ways to access this, including the OLS web API, SPARQL queries, or downloading the OWL and using an ontology API.

Comment as an answer, as requested.
The Incidence and Prevalance Database (IPD - subscription required).

license: Creative Commons - Attribution Share-Alike - https://creativecommons.org/licenses/by-sa/4.0/


text:
Why do reverse radtags have different start points in radtag sequencing?
Why is it that when doing paired end Radseq that R1 and R2 have different start points? For instance, while doing Radseq you do a restriction digestion of your DNA using enzyme X- shouldn't the sequencing reaction take place from either end of where the enzyme cut? If it doesn't- but happens randomly where do your multiplex sequences go?
For instance- if you are sequencing the sequence on an illumina machine with the "<>" representing the cut site-
Where does the sequencing reaction start from in each direction?
Where does the multiplexing sequence go? - does it go on each?
5-GGCTAGC<CCTGCA*GG>GCGCTT-3
 3-CCGATCG<GG*ACGTCC>CGCGAA-5

Hi @jawz_34, it appears you are asking several questions, and some of the questions are not relevant to the question title. It will be much easier for others to answer your post if you structure your post such that all of the post is remains relevant to the post title.

Assume that after digestion you're left with the following sequence:
5' AGCTAGCTAGCTAGCTGACGTCGTAC 3'
3' GTACGACGTCAGCTAGCTAGCTAGCT 5'

Then one read would start with AGCTAG and the other would start with TCGATC.
The barcodes (what are used for multiplexing) are outside of this and part of the illumina adapters. They're sequences in completely different reads. This is nicely illustrated below:

Index 1 and 2 don't actually need to be present. Even if they are, whether the sequencer actually sequences them is a user setting.
For those curious, the order of the reads as they're sequenced are:

read #1
Index 1 (optional)
Index 2 (optional)
read #2 (optional)

If you ever get more than 2 reads from a sequencing provider, then it's likely that they included the index sequence (I have no clue why, but some sequencing centers like to do this...even though it confuses the heck out of their users).


text:
Standard Cutoff for Moderated T-statistics
I'm looking at some microarray data. For the first time I've calculated a moderated T statistic from limma.
Is there any standard practice for where to cut off that value? For log2 fold change I usually cut off at +/- 1.5 and adjust accordingly; for adjusted p value I cut off at p<0.01 and adjust accordingly.
These moderated T statistics seem a little higher than log2 fold change and I'm not sure what is generally accepted as significant.
Thanks!

If you're filtering by the (adjusted) p-value, why do you want to also filter by moderated T-statistic? The p-value basically comes from the moderated T-statistic, after all.

I would like to have values that are real but are also above a particular magnitude. The P value gives me an idea of whether the finding is real, but I also need a way to cut-off on magnitude. I'd usually do this by filtering l|og2Fold|>1.5, but I was hoping to have a similar moderated T-statistic cutoff. Or am I misinterpreting the moderated T statistic?

You're misinterpreting the moderated T-statistic, it's basically the fold-change divided by its variance. The p-value comes directly from that, so if you filter by moderated fold-change you're just setting an unknown (unless you go through the trouble of figuring it out) p-value threshold. Instead, do exactly as you've been doing and set appropriate (adjusted) p-value and fold-change thresholds.


text:
Can exons be located outside of the coding sequence?
I have a gff file like this (I edited the name):
scaffold_x  source  exon    2987526 2987805 .   -   .   name "foobar";transcriptId 68892
scaffold_x  source  CDS     2987526 2987805 .   -   1   name "foobar";proteinId 68892;exonNumber 5
scaffold_x  source  exon    2988610 2988745 .   -   .   name "foobar";transcriptId 68892
scaffold_x  source  CDS     2988610 2988745 .   -   2   name "foobar";proteinId 68892;exonNumber 4
scaffold_x  source  exon    2989157 2989645 .   -   .   name "foobar";transcriptId 68892
scaffold_x  source  CDS     2989157 2989645 .   -   2   name "foobar";proteinId 68892;exonNumber 3
scaffold_x  source  exon    2990903 2991541 .   -   .   name "foobar";transcriptId 68892
scaffold_x  source  CDS     2990903 2991541 .   -   2   name "foobar";proteinId 68892;exonNumber 2
scaffold_x  source  stop_codon  2990903 2990905 .   -   0   name "foobar"
scaffold_x  source  exon    2991587 2992430 .   -   .   name "foobar";transcriptId 68892
scaffold_x  source  CDS     2991587 2992430 .   -   0   name "foobar";proteinId 68892;exonNumber 1
scaffold_x  source  start_codon 2992428 2992430 .   -   0   name "foobar"

So the start_codon is 2992428    2992430 and the stop_codon is 2990903   2990905; I would therefore expect all the exons within this range, however, e.g. exon 5 is found in 2987526   2987805 which is clearly outside this range.
Is there something wrong with the gff file itself or is there a plausible explanation for that?
EDIT:
In the corresponding protein fasta I do find an amino acid sequence associated with this prtein ID and name:
>68892|foobar
MHTGD...

Thanks for editing, but that doesn't tell us anything. We need to know the accession of the protein and/or gene you are looking at. What is 68892? Is that the gene ID? Also, why would there not be a protein sequence associated with it? It is only the UTRs that are not translated, the rest of the exons will be.

@terdon See my comment below your answer. And thanks for the hint; there are indeed exons in the "correct" region which can then be translated to a protein. Feel free to expand your answer and then I am happy to accept it.

I don't get what is confusing you, to be honest. I'd be happy to expand, but I don't know in what direction. It might be simpler if you come into chat for a second and ping me (@terdon) there so we can figure out what is confusing you.

Yes, of course. Exons are not limited to the protein coding regions. Many UTRs are in exons. In fact, you even have various cases of UTRs being multiple exons, and being spliced.
What is strange in your file is not so much that you have exons beyond the stop codon, but that you also have them marked as CDS (coding sequence). That isn't possible, no. While there can indeed be exons in the UTRs, those are not coding and shouldn't be called CDS. That's probably just a minor nomenclature issue though and what you are looking at is a spliced UTR.
So, the exons (or the parts of the exons) that fall between the start and stop codons will be translated, while those that fall outside that regions will not be (UTRs).
For example, this is the exonic structure shown for ENST00000617185, one of the transcripts of the human P53 gene:

The boxes are exons and the lines are introns. The colored boxes are protein coding exons, while the empty ones are UTR exons. Note how they are also spliced and that this particular transcript has 3 non-coding exons on the 5' end and one on the 3' end.

@Cleb no, UTRs are not translated. But I can't really help if you don't give us the gene you are looking at. Please edit your question and tell us where the data is coming from so we can also compare your gff to the predicted protein.

Thanks already for your answer (+ 1). So if I understand correctly, these regions should not encode for a protein, correct!? Funny thing is, however, if I look at the corresponding protein fasta, there is an amino acid sequence associated with it (see my edit). But, as it is not a coding sequence, this should not b the case then?

Unfortunately, it is not publicly available, when I google for the gene name I don't get any hits. But for now it is already good to know that there is something fishy(?) going on so that I can contact the authors; just wanted to make sure that I don't miss something obvious.

@Cleb there is nothing fishy going on at all. All (or, at least, the vast majority) of protein coding transcripts have untranslated regions (UTRs) and this is what you are showing. That doesn't mean that all of the transcript is not translated, only that the UTR regions are not translated. And what do you mean it isn't publicly available? Is this a gene prediction? All genes are publicly available.

Yes, thanks for clarifying, I am not very familiar with this topic, so I probably ask a couple of silly things :) Well, I got the file from a website where one can only download files if one is registered and as far as I understood one is not allowed to distribute the information without permission.

I'm wondering if the CDS outside of the start-stop range is to do with different isoforms? The exons have a transcript id, and the CDS entryies have a protein_id, but the start and stop lines only have gene ids. Could be that the CDS entries and the start/stop_codon entires are from different transcript isoforms of the same gene.

Those are the untranslated regions (UTRs). All mRNAs have a 5' UTR and a 3' UTR. These give the ribosome something to grab onto and often contain important regulatory sites such as miRNA target sites.

I've often wondered idly about that. Are you 100% sure that all mRNAs will have two UTRs? Or, at least, all mammalian mRNAs, to keep things simple. The vast majority do, of course, but this is biology, the land of shades of gray and exceptions. Are there absolutely no mRNAs without UTRs?

Not all have them identified (yet) but functionally I think that UTRs are necessary. I'm no expert on translation, but even ignoring the regulatory aspect, I think they're necessary just to allow the ribosome to attach on.

Yes, that's what I'm wondering about too. I'm not a translation expert either, but I don't remember any such function. That, however, in no way implies that there isn't one, which is why I thought I'd ask. I'm pretty sure I've seen transcripts with no annotated UTRs, but as you say, that might just be a fault in the annotation rather than an actual biological phenomenon.

Thanks already for your answer (+ 1). So if I understand correctly, these regions should not encode for a protein as it is a untranslated region, correct!? However, in the corresponding protein fasta file, I find an amino sequence associated with the name and protein ID (see my edit). How can this be?

Yes, the coding region encodes a protein. At either end of the coding region is a UTR.

Although the explanation that UTR (untranscribed) regions can consist of multiple exons covers most situations, I think it is good to mention stop codon readthrough.
This would result in a CDS (coding sequence) part of your mRNA (messenger RNA) after a stop codon.
However, following the gff3 format specification, this should probably be encoded differently, using SO (sequence ontology) terms like SO:0000697 (gene with stop codon read through) and  SO:0000883 (stop codon readthrough).

all_licenses: Creative Commons - Attribution Share-Alike - https://creativecommons.org/licenses/by-sa/4.0/


text:
What is the actual cause of excessive zeroes in single cell RNA-seq data? Is it PCR?
First, sorry if I am missing something basic - I am a programmer recently turned bioinformatician so I still don't know a lot of stuff. This is a cross post with a Biostars question hope that's not bad form.

While it is obvious that scRNA-seq data contain lots of zeroes, I couldn't find any detailed explanation of why they occur except for short notices along the lines of "substantial technical and biological noise". For the following text, let's assume we are looking at a single gene that is not differentially expressed across cells.
If zeroes were caused solely by low capture efficiency and sequencing depth, all observed zeroes should be explained by low mean expression across cells. This however does not seem to be the case as the distribution of gene counts across cells often has more zeroes than would be expected from a negative binomial model. For Example the ZIFA paper explicitly uses a zero-inflated negative binomial distribution to model scRNA-seq data. Modelling scRNA-seq as zero-inflated negative binomial seems widespread throughout the literature.
However assuming negative binomial distribution for the original counts (as measured in bulk RNA-seq) and assuming that every RNA fragment of the same gene from every cell has approximately the same (low) chance of being captured and sequenced, the distribution across single cells should still be negative binomial (see this question  for related math).
So the only remaining possible cause is that inflated zero counts are caused by PCR. Only non-zero counts (after capture) are amplified and then sequenced, shifting the mean of the observed gene counts away from zero while the pre-PCR zero counts stay zero. Indeed some quick simulations show that such a procedure could occasionally generate zero-inflated negative binomial distributions. This would suggest that excessive zeroes should not be present when UMIs are used - I checked one scRNA-seq dataset with UMIs and it seems to be fit well by plain negative binomial.
Is my reasoning correct? Thanks for any pointers.
The question
How can we distinguish between true zero and dropout-zero counts in single-cell RNA-seq? is related, but provides no clues to my present inquiry.

It may be necessary to distinguish between methods that use unique molecular identifiers (UMIs), such as 10X's Chromium, Drop-seq, etc, and non-UMI methods, such as SMRT-seq.  At least for UMI-based methods, the alternative perspective, that there is no significant zero-inflation in scRNA-seq, is also advocated in the single-cell research community.  The argument is straight-forward: the empirical mean expression vs. dropout rate curve matches the theoretically predicted one, given the current levels of capture efficiency.
Examples
Svensson Blog
A couple of blog posts from Valentine Svensson argue this point rather pedagogically, and include citations from across the literature:
Droplet scRNA-seq is not zero-inflated
Count-depth variation makes Poisson scRNA-seq data negative binomial
bayNorm
There is a more extensive preprint by Tang, Shahrezaei, et al. (BioRxiv, 2018) that claims to show a binomial model is sufficient to account for the observed dropout noise.  Here is a snippet of a relevant conclusion:

Importantly, as bayNorm recovered dropout rates successfully in both UMI-based and non-UMI protocols without the need for specific assumptions, we conclude that invoking zero-inflation models is not required to describe scRNA-seq data.  Consistent with this, the differences in mean expression levels of lowly expressed genes observed between bulk and scRNA-seq data, which were suggested to be indicative of zero-inflation, were recovered by our simulated data using the binomial model only.

Multinomial Modeling
There is also a very clearly written preprint by Townes, Irizarry, et al. (BioRxiv, 2019) where the authors consider scRNA-seq as a proper compositional sampling (i.e., multinomial process) and they come to a similar conclusion, though specifically for UMI-based methods. From the paper:

The multinomial model makes two predictions which we verified using negative control data. First, the fraction of zeros in a sample (cell or droplet) is inversely related to the total number of UMIs in that sample. Second, the probability of an endogenous gene or ERCC spike-in having zero counts is a decreasing function of its mean expression (equations provided in Methods). Both of these predictions were validated by the negative control data (Figure 1). In particular, the empirical probability of a gene being zero across droplets was well calibrated to the theoretical prediction based on the multinomial model.  This also demonstrates that UMI counts are not zero inflated.

Furthermore, by comparing raw read counts (prior to UMI-based deduplication) and UMI counts, they conclude that PCR is indeed the cause of zero-inflation:

The results suggest that while read counts appear zero-inflated and multimodal, UMI counts follow a discrete distribution with no zero inflation (Figure S1). The apparent zero inflation in read counts is a result of PCR duplicates.

I highly recommend giving this a read, especially because it nicely situates other common generative models (e.g., binomial, Poisson) as valid simplifying assumptions of the multinomial model.
It should be noted that this same group previously published a work (Hicks, Irizarry, et al. 2018), mostly focused on non-UMI-based datasets (SMRT-seq), where they showed evidence that, relative to bulk RNA-seq, there was significant zero-inflation.

One more preprint arguing that zero inflation is not a good model: https://www.biorxiv.org/content/10.1101/477794v1.abstract

The Biostars thread turned out helpful. The most interesting possible cause, not mentioned in the Ian Subery's answer, is that due to bursty nature of transcription, the true distribution of transcript counts across cells can be bimodal with a peak at zero even assuming a simple model of transcription such as the random telegraph model. See for example Dattani & Barahona 2017 for a more detailed discussion.
However, the random telegraph model predicts that the even with peak at zero, there should usually be non-negligible probability of having small, but nonzero counts, which is not always the case in zero-inflated scRNA-seq data: see for example the Isl1 gene in Retina dataset (click Explore -> search for "Isl1")

I know of no references for this, but in general, I would say that your reasoning is sound. I would just add that in contrast to what I suspect you have simulated, not all transcripts are equally likely to be captured and amplified. We don't really understand what the determinants of this are, but for example, GC content is definitely related.

Thanks for the note. I understand there are betweeen-gene differences due to GC content etc. but I am only looking at one gene at a time so it shoud not matter. Do you believe that the sequencing results could be noticeably influenced by within-gene differences - e.g. because a specific fragment of a long transcript is much better captured/amplified than others letting us observe large non-zero counts only in cells where this fragment was captured and zeroes otherwise?


text:
Why are inversions defined as the reverse complement and not just the reverse of the reference?
I can’t quite understand the way inversions are defined. In particular I expect an inversion to be only a reversed version of the reference and not its reverse complement. Most sources use diagrams such as below which totally obscure the fact that it is also complemented.

Could someone provide an explanation on what causes an inversion and why it is defined the way it is?

If an inversion were just reversed then there would be 3' -> 3' bonds and 5' -> 5' bonds. That's why inversions are reverse complemented, you then maintain the normal 5'->3' direction. This is nicely covered over on our biology sister site, which also includes this image:

Thanks for the link. The image you posted made it clear for me. The key is to look at both strands at the same time.


text:
How to get UCSC Genome Browser to show the actual sequence?
Is it possible to make UCSC Genome Browser to show the sequence (nucleotide chain) of the current region?
I fiddled with the various viewing options and couldn't figure it out.

If you zoom in at one specific region, it shows the nucleotide chain on the third line (below the chromosome number). It can't be there by default as the size of default region is more than 100,000bp. You can use the scale on the top of the page to see the size of the region.

In case it’s of any benefit to you, I’ll often check IGV when looking at tracks with that kind of zoomage. Not that it zooms any better than the UCSC browser (although the bases are colour coded!) but you can see your aligned reads individually if you want, along with directional information etc etc. Loads of stuff.

From the top bar:
View -> DNA


text:
Genome assembly from error-prone reads
I understand how to assemble genome from error-free reads. I implemented like this:

Construct directed overlap graph with reads as vertices and edges as
  maximum overlap between two vertices. Edges represent the length of
  overlapping maximum substring of one read suffix with other read
  prefix. Find the path that visits all vertices only once in overlap
  graph. With small graph we can use greedy approach : Starting from
  first pattern always choose egde of maximum weight and store shortest
  common superstring spelled by this path in resulting string. Resulting
  string spelled by path is shortest common superstring of all reads,
  that is assembled genome.

For example
AACTAG
  CTAGAT

Gives the AACTAGAT as a result
Now let's say that two error-prone reads have at most one error (mismatch):
AACTAG
  CTCGAT

As you can see second read has one substitution in third character from 'A' to 'C'.
So how assembled genome should look like for the error-prone reads? Do I take A or C to resulting string? Or I'm missing something?
Note: I used short reads for example because of simplicity. I work with much larger number of reads which are 100 nucleotides long and have at most one substitution of character. So that's a big overlap graph. I know I can use De Bruijn graph, but can I use some greedy approach (like in error-free example)?

It's common for assemblers to have a coverage threshold to make sure that there are enough reads across the same region to properly assemble a contig.
You can't work out an ideal solution from a situation where the consensus sequence is ambiguous. If your overlap were representative (e.g. a heterozygous variant in a diploid genome, based on a few tens of reads of 100bp, rather than 6bp), then the resulting assembly would have a bubble in the assembly graph at the variant point.
        A
       / \
A-A-C-T   G-A-T
       \ /
        C

A conservative assembler might produce two assembled contigs corresponding to these two paths, AACTAGAT and AACTCGAT. A greedy assembler would just pick one path and use that. Subsequent assembly steps might attempt to search for bubbles and annotate the assembled contigs with variant locations.

For all mainstream short-read overlappers, if there is a mismatch, the overlap will be missed entirely. You don't get a bubble. Usually you need at least 4 reads to see a bubble and the bubble is not single-based in string/overlap graph.

Yes, for the particular example that was given, a bubble wouldn't be formed. I was assuming that the reads were representative of something with a better coverage (e.g. overlap of 50 long sequences, with half having a variant, and half not)

For short reads, the typical and the most widely used solution is to correct away sequencing errors before assembly. You can correct errors with k-mer spectrum, a trie or multi-alignment. There are many papers on this topic.
Error correction alone won't fix all sequencing errors. Remaining errors may lead to bubbles and tips in the overlap/de Bruijn graph. You can identify them by inspecting the local topology of assembly graph.

So what would be solution from error-prone sample above?

No solution to your example. You can only identify errors with sufficient coverage.


text:
Which identifier do I have to use when I want to add a ##Fasta section to a gff?
I have a gff file like this
scaffold_1  source  mRNA    2987526 2992430 .   -   .   ID=protein_id_68892;Name="foobar";locus_tag=my_organism_68892;translation=MHTGDALEGSTGNVSILV...
scaffold_1  source  exon    2987526 2987805 .   -   .   name "foobar"; transcriptId 68892;locus_tag=my_organism_68892
scaffold_1  source  CDS     2987526 2987805 .   -   1   name "foobar"; proteinId 68892; exonNumber 5;locus_tag=my_organism_68892;Parent=protein_id_68892
scaffold_1  source  exon    2988610 2988745 .   -   .   name "foobar"; transcriptId 68892;locus_tag=my_organism_68892
scaffold_1  source  CDS     2988610 2988745 .   -   2   name "foobar"; proteinId 68892; exonNumber 4;locus_tag=my_organism_68892;Parent=protein_id_68892
scaffold_1  source  exon    2989157 2989645 .   -   .   name "foobar"; transcriptId 68892;locus_tag=my_organism_68892
scaffold_1  source  CDS     2989157 2989645 .   -   2   name "foobar"; proteinId 68892; exonNumber 3;locus_tag=my_organism_68892;Parent=protein_id_68892
scaffold_1  source  exon    2990903 2991541 .   -   .   name "foobar"; transcriptId 68892;locus_tag=my_organism_68892
scaffold_1  source  CDS     2990903 2991541 .   -   2   name "foobar"; proteinId 68892; exonNumber 2;locus_tag=my_organism_68892;Parent=protein_id_68892
scaffold_1  source  exon    2991587 2992430 .   -   .   name "foobar"; transcriptId 68892;locus_tag=my_organism_68892
scaffold_1  source  CDS     2991587 2992430 .   -   0   name "foobar"; proteinId 68892; exonNumber 1;locus_tag=my_organism_68892;Parent=protein_id_68892

The mRNAentry was not present at first and I added it in. I also have the protein fast file available and would now like to add the amino acid sequences to the gff. I originally added them in column 9 to the mRNAentry which is apparently incorrect. Here it is written that the sequences should go at the end of the file in the ##Fastasection. However, I struggle to figure out which identifier I should use to link it to my mRNA entries.
Do I have to use
##Fasta
>protein_id_68892
MHTGDALEGSTGNVSILV...

or does it work with the locus_tag
##Fasta
>my_organism_68892
MHTGDALEGSTGNVSILV...

or is it something completely else? Where and how do can I add the translation information?
In the link above (and also in the comments it is suggested) the identifier in the first column is used, in my case that would be scaffold_1. However, I do have a lot of entries with scaffold_1in my file so a link between scaffold_1and an amino acid sequence would not be unique as I would have
##Fasta
>scaffold_1
MHTGDALEGSTGNVSILV...
>scaffold_1
GDAHTEGSTGNVSILV...

. In the end, I want to create a genbank file (which does not exist yet) where I can link a CDS to an amino acid sequence.

Normally, the ##Fasta section contains the nucleotide sequences corresponding to the seqid column. Most GFF parsers can generate amino acid sequences on the fly from this information

@heathobrien: Thanks for the info. But the problem would remain: Which identifier should be used there?

in your example, scaffold_1

@heathobrien: but in the file I see scaffold_1 for a lot of proteins.  Then it would still unclear which sequence to add. Also, in the link in my question is written This notation indicates that the annotation portion of the file is at an end and that the remainder of the file contains one or more sequences (nucleotide or protein) in FASTA format. So adding proteins should be ok. Any ideas?

It's normal for multiple genes to be coded on the same DNA sequence. The GFF tells you which portion of the sequence codes for each gene. Ultimately, you can add whatever you want. It really depends what you're planning to do with the file.

Unless I'm mistaken, the only variant of GFF that explicitly permits Fasta entries is GFF3, and the data you've shown above is definitely not valid GFF3. The mRNA entry seems to correctly use the = symbol to separate keys from values in each attribute, but the remaining columns all use spaces and quotation marks (GTF style) instead. This is likely to cause all kinds of problems with software packages that take GFF data as input.
But to more directly address your question, there is no "right" way to connect feature entries to sequence entries in a GFF3 file. There are some conventions, but these are all very specific to a small set of tools and community of practice, and it is rare that these conventions are objectively enforced. Some options include:

Placing genomic nucleotide sequences in the ##FASTA section, and then using feature coordinates to pull out mature mRNA sequences and translate them in silico to derive the corresponding peptide sequence. In this case, the Fasta entries would be genomic sequence identifiers that correspond to the first field/column in each feature entry.
Placing mature mRNA sequences in the ##FASTA section, and then directly translating this sequence to derive the corresponding peptide sequence. In this case, you'll want the sequence IDs in the ##FASTA section to match the ID attribute of the corresponding mRNA feature.
Placing the peptide sequences directly in the ##FASTA section. In this case, you would either want to have protein features declared in the GFF3 file (uncommon) or have some sort of protein_id attribute associated with each mRNA feature, in which case you'd want the sequence IDs in the ##FASTA section to match this protein_id attribute.

There are probably other ways of doing it, too. Which way is "right"? The GFF3 specification has nothing to say about it. It depends on the tools you are using to do the subsequent processing and analysis. If you're lucky, those tools might even document what they expect in terms of input.

Let me just say, as a post-script, that I think placing sequences and feature annotations in the same file is a horrible idea. I've never liked this about GFF3.

To answer your specific question, you should use the identifier from the first column of the GFF in the ##Fasta section. As you note, this won't work if you are trying to add multiple protein sequence for each seqid. Try using the scaffold DNA sequences in this section and add amino acid translations as an attribute of the annotations.

Thanks for your reply. add amino acid translations as an attribute of the annotations Which annotations? Where exactly would you add the translations?


text:
Counting number of possible alignments between two DNA sequences using python
I was looking BioPython to calculate the number of possible alignments between two sequences. Let's say there are two input sequences with the length of m and n. The program will count the possible number of alignments based on the length of two sequences.
Is it possible to do that with BioPython? If not, how can I do that?

Which kind of sequences DNA, RNA or proteins? How have you tried to do so? What do you expect to learn from it?

And the third question? Why do you need to know? (By the way, if you don't try why do you expect random people on the internet to help you?, Not that someone can't but it is harder to provide help if one doesn't see research effort)

Seems like a homework assignment where they are asking for the theoretical number of alignments.  even down to the typical n and m nomenclature.

You should define what you mean by "possible alignment". Biopython will only help you use tools that chose among the possible alignments, based on certain criteria.

A closed-form solution is offered in An exact formula for the number of alignments
between two DNA sequences by Torres, Cabada, and Nieto:
$$f(m,n)=\sum_{k=0}^{min(m,n)}2^{k}\binom{m}{k}\binom{n}{k}$$
If this solution seems reasonable, you could calculate this without BioPython, but with a simple xrange loop, the power operator, and scipy.special.binom:
#!/usr/bin/env python

import sys
import scipy.special

def alignments(m, n):
    if m < 0 or n < 0:
        raise ValueError('m and n should be non-negative')
    s = 0
    b = min(m, n) + 1
    for k in xrange(0, b):
        s += (2 ** k) * scipy.special.binom(m, k) * scipy.special.binom(n, k)
    return s

def main():
    m = 4
    n = 2
    a = alignments(m, n)
    assert (m == 4 and n == 2 and a == 41), "something went wrong!"
    sys.stdout.write("alignments(%d, %d) = %d\n" % (m, n, a))

if __name__ == "__main__":
    main()

I don't know what your m and n are, but keep in mind that this solution grows quickly. Also, I did some digging and scipy approximates the calculation of the binomial coefficient, which is probably going to be an issue for large m and n.
If you're working with large values, you may want to look into using the Python long type and consider the accuracy of the answer. The Stack Overflow link offers alternatives to scipy that generate an exact calculation of the binomial coefficient, at the cost of a longer runtime.


text:
Chimera Alignments
I have a structure with two subunits. I am trying to show movement of the C-terminal subunit upon ligand binding by superposition with another structure from the same strain in the apo form.
I want to superpose the N-terminal subunits (A), note the root mean square of the atomic deviations (RMSD), then apply the resulting matrix to the superposition of both subunits of the molecule with both subunits of the other (B), then I want to superpose on the C-terminal subunits (C), and then measure the RMSD between (C) and the C-terminal subunits from (B).
Do you know how can I do this in Chimera? I know how to superpose molecules and calculate the RMSD in Chimera, but I do not know how to apply the resulting matrix from (A) to (B). I tried to save the aligned structures from (A) as a PDB file to check and see if there was a matrix in the file, but an error message came up saying that I needed to save the file as $name or $number. I tried to save the file with $ in front of the name, but it didn't work.
How can I do this correctly?

Answer from @terdon converted from comment:
The $ are often used to signify that something is a variable. So I guess they mean you need to use a name that you've already used. I mean, the same name as used for something else in your pipeline, so that then the program will see the file with the same name and know to apply it. Sorry, I don't have any experience with chimera, so I have no idea if what I'm suggesting is reasonable. It's just what I get from the $ you mentioned.
Answer from @marcin converted from comment:
this question was answered on chimera-users:
http://plato.cgl.ucsf.edu/pipermail/chimera-users/2017-October/014045.html
Answer from @joe-healey converted from comment:
You are probably looking for the “Morph Conformations” function within Chimera. It’ll move and superimpose 2 structures on to one another.


text:
Creating a new reference genome B by changing genome A with mismatches from BAM file without long reads
This question is slightly related to this one:
Improve a reference genome with sequencing data
Only in my case, I have a starting genome reference A, of a short phage genome, and reads from a related genome B, for which I haven't been able to find a reference.
I took the Illumina 2x75bp paired end reads from genome B and ran Velvet on it (I ran Ray, but it didn't complete). It gives me a fragmented collection of contigs. The contigs align 99-98% to reference genome A. I don't have long reads.
My question is: with the short reads and the contigs that I have, can I transform the reference genome A into a reference genome B?
PS: no, I don't have long reads, and there is no budget for them. The budget for long reads is zero and will remain zero. Even if we could sequence long reads for free, we would still not do it because of the lab work.
Tools that I have tried so far:
ococo:
$HOME/ococo/ococo -i contigs.bwamem.bam -c 1 -f $HOME/Genomes/T2_phage/starting.fasta -F consensus.fasta --vcf-cons - -V variants.vcf

It produced a consensus.fasta, and checking how much it changes from the reference, it changes. There are one or more changes in about 33% of the 100bp windows:
diff -y <(seqtk seq -l 100 $HOME/Genomes/T2_phage/starting.fasta) <(seqtk seq -l 100 consensus.fasta) | grep '\ \|\ ' | wc -l
1272
diff -y <(seqtk seq -l 100 $HOME/Genomes/T2_phage/starting.fasta) <(seqtk seq -l 100 consensus.fasta) | grep -v '\ \|\ ' | wc -l
419

What is the insert size distrobution on the 2x75 library? Based on the alignments what is the coverage of the reference A genome?  How does the N50 of the velvet assembly compare with the reference A genome size?

The size distribution is 250bp. The alignment coverage is about 200-250x. The genome size of the starting genome is about 157kb, the two longest contigs were about 2-2.5kb.

Why aren't you using the solution to the question you linked to? It seems very similar to your situation. (If you describe the difference explicitly, it might help you get better answers on this thread.)

Even though ococo does not seem to handle big indels well, it's a good enough tool for a first try.

You should be able to use something like Pilon to convert from one genome to another, assuming the changes are all fairly small local changes (but that's almost never the case).
I find your resistance to long-read data interesting. It's important to consider the cost of bioinformatics as well. Sequencing a phage genome with long reads will take less than an hour, and probably produce a few reads that cover the entire genome. With existing short-read data, those full-length reads should be able to be fixed up within a day.
I spent about half a year trying to fix a eukaryotic genome assembly that was assembled only from short reads (using Pilon, and L-RNA-scaffolder), and I was only able to reduce the number of contigs down by about half (and even at that, I wasn't convinced that contigs were being properly joined). With long reads and another couple of weeks of computer work, the contigs dropped down to about a tenth of their original number.
If you can use a pipette and have about half an hour spare, you can do rapid long-read sequencing on the MinION.


text:
How do I export my biom object to a file in biom format?
I'm using the python API for the biom-format package. I read in two files in biom format and merged them, like so:
In [1]: import biom

In [2]: a = biom.load_table('./table-a.biom')

In [3]: b = biom.load_table('./table-b.biom')

In [4]: a
Out[4]: 51116 x 46 <class 'biom.table.Table'> with 185077 nonzero entries (7% dense)

In [5]: b
Out[5]: 37479 x 26 <class 'biom.table.Table'> with 100124 nonzero entries (10% dense)

In [6]: c = a.merge(b)

In [7]: c
Out[7]: 67022 x 47 <class 'biom.table.Table'> with 254555 nonzero entries (8% dense)

How do I now export this biom object (c) as a biom file? There are a few methods available, eg, .to_tsv(), but I don't see a way to save this object in biom format.

Are you using the biomformat package of bioconductor through Python? Or the biom-format package? Anyway in R the relevant function is write_biom so it might be the same in the python package

You can use the functions to_hdf5 and to_json from the biom-format package:
to_hdf5
with biom.util.biom_open('output.biom', 'w') as f:
    merged_table.to_hdf5(f, "example")

to_json
with open('output.txt', 'w') as f:
    f.write(merged_table.to_json("example"))

There is a function in PiCRUST you can use
from picrust.util import write_biom_table

write_biom_table(c,'table-c.biom')

Note, however, that depending if you want json or hdf5 format and if your metadata are complex, you may need to use additional arguments. These are the arguments of the function:
def write_biom_table(biom_table, biom_table_fp, compress=True,
                 write_hdf5=HAVE_H5PY, format_fs=None):
"""Writes a BIOM table to the specified filepath

Parameters
----------
biom_table : biom.Table
    The table object to write out
biom_table_fp : str
    The path to the output file
compress : bool, optional
    Defaults to ``True``. If True, built-in compression on the output HDF5
    file will be enabled. This option is only relevant if ``write_hdf5`` is
    ``True``.
write_hdf5 : bool, optional
    Defaults to ``True`` if H5PY is installed and to ``False`` if H5PY is
    not installed. If ``True`` the output biom table will be written as an
    HDF5 binary file, otherwise it will be a JSON string.
format_fs : dict, optional
    Formatting functions to be passed to `Table.to_hdf5`

Notes
-----
This code was adapted from QIIME 1.9
"""

I have some problems to perform a simple operation and write the output of biom data generated with picrust, it might be helpful as an example of how to use the arguments.


text:
How to map PDB chains to Uniprot IDs using API services
I have a lot of PDB IDs and I need to get uniprot fasta sequences of these PDB IDs special chains by API services. For example, imagine that I need to get fasta sequence of '1kf6' 'A' chain. The uniprot entry (Accession) for this '1kf6' chain is 'P00363'. If I find a way to get this entry, I can simply retrieve its fasta sequence using bioservices package:
from bioservices import *
u = UniProt()
sequence = u.retrieve("P00363","fasta")
print(sequence)

The main problem is that I don't know how can I find this uniprot entry("P00363" here) for each chains in PDB IDs(here I imagine that I don't know this "P00363" is the uniprot entry for '1kf6' 'A' chain and I'm trying to find it using API services). I tried to use ID mapping. In ID mapping cross-references, the chain name for each entry is specified and if I can add a cross-reference column to u.search('1kf6', columns="id,entry name,length, genes,pdb"), I can simply find the entry related to my favourite chain. The following code returns some entries each related to one of '1kf6' chains. The problem is that I don't know exactly which entry is related to 'A' chain that here is my favorite chain. In the following code, id searches for entry and other parameters search for sth else related to them:
from bioservices import *
u = UniProt()
res = u.search("1kf6",columns="id,entry name,length, genes")
print(res)

Each one of the following words can be used and added to columns part:
'id', 'entry name', 'genes', 'genes(PREFERRED)', 'genes(ALTERNATIVE)', 'genes(OLN)', 'genes(ORF)', 'organism', 'organism-id', 'protein names', 'proteome', 'lineage(ALL)', 'lineage-id', 'virus hosts', 'fragement', 'sequence', 'length', 'mass', 'encodedon', 'comment(ALTERNATIVE PRODUCTS)', 'comment(ERRONEOUS GENE MODEL PREDICTION)', 'comment(ERRONEOUS INITIATION)', 'comment(ERRONEOUS TERMINATION)', 'comment(ERRONEOUS TRANSLATION)', 'comment(FRAMESHIFT)', 'comment(MASS SPECTROMETRY)', 'comment(POLYMORPHISM)', 'comment(RNA EDITING)', 'comment(SEQUENCE CAUTION)', 'feature(ALTERNATIVE SEQUENCE)', 'feature(NATURAL VARIANT)', 'feature(NON ADJACENT RESIDUES)', 'feature(NON STANDARD RESIDUE)', 'feature(NON TERMINAL RESIDUE)', 'feature(SEQUENCE CONFLICT)', 'feature(SEQUENCE UNCERTAINTY)', 'version(sequence)', 'domains', 'domain', 'comment(DOMAIN)', 'comment(SIMILARITY)', 'feature(COILED COIL)', 'feature(COMPOSITIONAL BIAS)', 'feature(DOMAIN EXTENT)', 'feature(MOTIF)', 'feature(REGION)', 'feature(REPEAT)', 'feature(ZINC FINGER)', 'ec', 'comment(ABSORPTION)', 'comment(CATALYTIC ACTIVITY)', 'comment(COFACTOR)', 'comment(ENZYME REGULATION)', 'comment(FUNCTION)', 'comment(KINETICS)', 'comment(PATHWAY)', 'comment(REDOX POTENTIAL)', 'comment(TEMPERATURE DEPENDENCE)', 'comment(PH DEPENDENCE)', 'feature(ACTIVE SITE)', 'feature(BINDING SITE)', 'feature(DNA BINDING)', 'feature(METAL BINDING)', 'feature(NP BIND)', 'feature(SITE)', 'go', 'go(biological process)', 'go(molecular function)', 'go(cellular component)', 'go-id', 'interpro', 'interactor', 'comment(SUBUNIT)', 'citation', 'citationmapping', 'created', 'last-modified', 'sequence-modified', 'version(entry)', '3d', 'feature(BETA STRAND)', 'feature(HELIX)', 'feature(TURN)', 'comment(SUBCELLULAR LOCATION)', 'feature(INTRAMEMBRANE)', 'feature(TOPOLOGICAL DOMAIN)', 'feature(TRANSMEMBRANE)', 'annotation score', 'score', 'features', 'comment(CAUTION)', 'comment(TISSUE SPECIFICITY)', 'comment(GENERAL)', 'keywords', 'context', 'existence', 'tools', 'reviewed', 'feature', 'families', 'subcellular locations', 'taxonomy', 'version', 'clusters', 'comments', 'database', 'keyword-id', 'pathway', 'score'

I used some of the above words in columns part but they couldn't retrieve chain name. Anyone know which one of the above parameters can give me the chain name for this PDB ID? or anyone has any other idea for finding uniprot entry(UniprotKB AC) related to each PDB ID chain by API services?

I found the answer and I thought it may be useful for others. we can use pypdb package for that purpose (it works in unix systems):
import pypdb

all_info = pypdb.get_all_info('1kf6')
print(all_info)

the output is as follows:
{'polymer': [{'@entityNr': '1', '@length': '602', '@type': 'protein', '@weight': '66057.6', 'chain': [{'@id': 'A'}, {'@id': 'M'}], 'Taxonomy': {'@name': 'Escherichia coli', '@id': '562'}, 'macroMolecule': {'@name': 'Fumarate reductase flavoprotein subunit', 'accession': {'@id': 'P00363'}}, 'polymerDescription': {'@description': 'FUMARATE REDUCTASE FLAVOPROTEIN'}, 'enzClass': {'@ec': '1.3.5.4'}}, {'@entityNr': '2', '@length': '243', '@type': 'protein', '@weight': '27021.9', 'chain': [{'@id': 'B'}, {'@id': 'N'}], 'Taxonomy': {'@name': 'Escherichia coli', '@id': '562'}, 'macroMolecule': {'@name': 'Fumarate reductase iron-sulfur subunit', 'accession': {'@id': 'P0AC47'}}, 'polymerDescription': {'@description': 'FUMARATE REDUCTASE IRON-SULFUR PROTEIN'}, 'enzClass': {'@ec': '1.3.5.1'}}, {'@entityNr': '3', '@length': '130', '@type': 'protein', '@weight': '14898.8', 'chain': [{'@id': 'C'}, {'@id': 'O'}], 'Taxonomy': {'@name': 'Escherichia coli', '@id': '562'}, 'macroMolecule': {'@name': 'Fumarate reductase subunit C', 'accession': {'@id': 'P0A8Q0'}}, 'polymerDescription': {'@description': 'FUMARATE REDUCTASE 15 KDA HYDROPHOBIC PROTEIN'}}, {'@entityNr': '4', '@length': '119', '@type': 'protein', '@weight': '13118.9', 'chain': [{'@id': 'D'}, {'@id': 'P'}], 'Taxonomy': {'@name': 'Escherichia coli', '@id': '562'}, 'macroMolecule': {'@name': 'Fumarate reductase subunit D', 'accession': {'@id': 'P0A8Q3'}}, 'polymerDescription': {'@description': 'FUMARATE REDUCTASE 13 KDA HYDROPHOBIC PROTEIN'}}], 'id': '1KF6'}

This output gives all the information related to all chains and we can simply find uniprotkb AC of each chain in this output. I also found another answer for this question. This answer can be found here and the following code:
import requests
from xml.etree.ElementTree import fromstring

pdb_id = '1kf6.A'
pdb_mapping_url = 'http://www.rcsb.org/pdb/rest/das/pdb_uniprot_mapping/alignment'
uniprot_url = 'http://www.uniprot.org/uniprot/{}.xml'

def get_uniprot_accession_id(response_xml):
    root = fromstring(response_xml)
    return next(
        el for el in root.getchildren()[0].getchildren()
        if el.attrib['dbSource'] == 'UniProt'
    ).attrib['dbAccessionId']

def get_uniprot_protein_name(uniport_id):
    uinprot_response = requests.get(
        uniprot_url.format(uniport_id)
    ).text
    return fromstring(uinprot_response).find('.//{http://uniprot.org/uniprot}recommendedName/{http://uniprot.org/uniprot}fullName'
    ).text

def map_pdb_to_uniprot(pdb_id):
    pdb_mapping_response = requests.get(
        pdb_mapping_url, params={'query': pdb_id}
    ).text
    uniprot_id = get_uniprot_accession_id(pdb_mapping_response)
    uniprot_name = get_uniprot_protein_name(uniprot_id)
    return {
        'pdb_id': pdb_id,
        'uniprot_id': uniprot_id,
        'uniprot_name': uniprot_name
    }

print (map_pdb_to_uniprot(pdb_id))

The output of this code is as follows:
{'pdb_id': '1kf6.A', 'uniprot_id': 'P00363', 'uniprot_name': 'Fumarate reductase flavoprotein subunit'}

An alternative solution involves the use of the SIFTS resource from the EBI. I wrote a python parser for their xml format, which will give you a mapping between a pdb and a uniprot entry at the single residue level.
$ curl --silent "ftp://ftp.ebi.ac.uk/pub/databases/msd/sifts/xml/1kf6.xml.gz" | gunzip | python parse_sifts.py
1kf6    A       M       0       P00363  M       1
1kf6    A       Q       1       P00363  Q       2
1kf6    A       T       2       P00363  T       3
1kf6    A       F       3       P00363  F       4
1kf6    A       Q       4       P00363  Q       5
1kf6    A       A       5       P00363  A       6
1kf6    A       D       6       P00363  D       7
1kf6    A       L       7       P00363  L       8
1kf6    A       A       8       P00363  A       9
1kf6    A       I       9       P00363  I       10
[...]

The returned fields are as follows:

pdb ID
chain
pdb residue
pdb position
uniprot ID
uniprot residue
uniprot position

In the above snippet, just put your PDB ID of choice to retrieve the data you need.

It seems to be good but when I entered the code ($ curl --silent "ftp://ftp.ebi.ac.uk/pub/databases/msd/sifts/xml/1kf6.xml.gz" |                    gunzip | python parse_sifts.py) in Cygwin, I got the following error: python: can't open file 'parse_sifts.py': [Errno 2] No such file or directory
(23) Failed writing body

Are you sure that the Python script was in the same directory where you executed the snippet?

When I try to open "ftp.ebi.ac.uk/pub/databases/msd/sifts/xml/1kf6.xml.gz" in my browser, I'll get this error: This site can’t be reached
The webpage at ftp://ftp.ebi.ac.uk/pub/databases/msd/sifts/xml/1kf6.xml.gz" might be temporarily down or it may have moved permanently to a new web address.
ERR_FTP_FAILED.
I think the problem is because of the server.

I just tried and it worked for me fine: ftp://ftp.ebi.ac.uk/pub/databases/msd/sifts/xml/1kf6.xml.gz

Now it just worked for me too...

While not as nice as using SIFT. You can use the uniprot sparql endpoint as well to get mappings.
PREFIX up: <http://purl.uniprot.org/core/>
SELECT
  ?pdbId ?chain ?pdbChain ?uniprot
WHERE
{
  # A space separated list of pairs of PDB identifiers and the chain code.
  VALUES(?pdbId ?pdbChain) { ('6VXC' 'A') ('1BG3' 'B') }

  # Make an IRI out of the pdbId
  BIND(iri(concat('http://rdf.wwpdb.org/pdb/', ?pdbId)) AS ?pdb)

  # Map to UniProt entries
  ?uniprot rdfs:seeAlso ?pdb .
  ?pdb up:database <http://purl.uniprot.org/database/PDB> ;
   up:chainSequenceMapping ?chainSm .
  ?chainSm up:chain ?chainsPlusRange .

  # Extract the list of chains from the text representation.
  BIND(STRBEFORE(?chainsPlusRange, '=') AS ?chain)

  # Filter those that match.
  FILTER(CONTAINS(?chain, ?pdbChain))
}

Check out this file on SIFT.
ftp.ebi.ac.uk/pub/databases/msd/sifts/csv/pdb_chain_uniprot.csv
It contains residue-wise mapping of PDB chains with uniprot ids.
# 2020/04/12 - 21:03 | PDB: 15.20 | UniProt: 2020.02
PDB,CHAIN,SP_PRIMARY,RES_BEG,RES_END,PDB_BEG,PDB_END,SP_BEG,SP_END
101m,A,P02185,1,154,0,153,1,154
102l,A,P00720,1,40,1,40,1,40
102l,A,P00720,42,165,41,None,41,164
.
.


text:
Mark BWA-SW split alignments in output for long reads
Is there a way to remove split alignments of single-end long reads from BWA-SW output? BWA-MEM has an option to include or flag split alignments in the result, but it seems BWA-SW method include them in the result unless the split occur in paired-end reads.

From BWA manual, it is advised to use BWA-MEM in your longs read case (unless you have a good reason not to):

BWA-MEM and BWA-SW share similar features such as long-read support and split alignment, but BWA-MEM, which is the latest, is generally recommended for high-quality queries as it is faster and more accurate.

I guess that BWA-SW returns several lines per transcript (based on local alignments). In that case, the issue has already been reported. The user asks:

How can I differentiate between a read that has multiple hits in the reference and a read that has a split alignment?

The author replies:

In principle, it is hard to tell. Bwa-mem makes an arbitrary choice.
  Mappings flagged with 0x100 are multiple hits. Mappings flagged with
  0x800 are from a chimeric alignment.

So you could actually filter by flag. But I am afraid the output could include something you do not necessarily want to (since there is no specific flag for the split reads).
Maybe a better solution: you could switch to some mapper that will report the proper split in the hit of the read map, such as minimap2 from the same author. It is suitable for long reads and comes with a splice -x feature.
Then, the laziest solution I see here would be to remove the split alignments exploiting the sam format and in particular the CIGAR string. Meaning, to remove the sam fields containing a N in their CIGAR (aka, skipped region from the reference). Something like this would work:
cat alignment.sam | awk '$6 !~ /N/ {print $0}'

Hope this can be somewhat helpful.


text:
Finding out proteins that have a similar domain to a specific toxin
I am trying to find out that which DNA sequences(in the below) have a similar domain to a specific toxin which is Clostridium botulinum
You can find the fragment of the DNA sequence that belong to that toxin gene below:

gi|27867582 (fragment of the known Clostridium botuninum toxin gene)
  GTGAATCAGCACCTGGACTTTCAGATGAAAAATTAAATTTAACTATCCAAAATGATGCTT
  ATATACCAAAATATGATTCTAATGGAACAAGTGATATAGAACAACATGATGTTAATGAAC
  TTAATGTATTTTTCTATTTAGATGCACAGAAAGTGCCCGAAGGTGAAAATAATGTCAATC
  TCACCTCTTCAATTGATACAGCATTATTAGAACAACCTAAAATATATACATTTTTTTCAT
  CAGAATTTATTAATAATGTCAATAAACCTGTGCAAGCAGC

I ran a blastx job over the above sequence which translates nucleotide into protein. Here's the result.
The result tells me that. The particular DNA sequence listed above, belongs to botulinum neurotoxin type E [Clostridium botulinum]
Besides, I am given two DNA sequences:
>Sequence A
TCTATCAAGTAGATTATTAAATACTACTGCACAAAATGATTCTTACGTTCCAAAATATGA
TTCTAATGGTACAAGTGAAATAAAGGAATATACTGTTGATAAACTAAATGTATTTTTCTA
TTTATATGCACAAAAAGCTCCTGAAGGTGAAAGTGCTATAAGTTTAACTTCTTCAGTTAA
TACAGCATTATTAGATGCATCTAAAGTTTATACGTTTTTTTCTTCAGATTTTATTAATAC

> Sequence B
TCCTGGCTCAGGACGAACGCTGGCGGCGTGTGCTTAACACATGCAAGTCGAGCGATGAAG
CTTCCTTCGGGAAGTGGATTAGCGGCGGACGGGTGAGTAACACGTGGGTAACCTGCCTCA
AAGTGGGGGATAGCCTTCCGAAAGGAAGATTAATACCGCATAACATAAGAGAATCGCATG
ATTTTCTTATCAAAGATTTATTGCTTTGAGATGGACCCGCGGCGCATTAGCTAGTTGGTA

Again, I translated them into proteins. The result are listed respectively(A and B):
https://blast.ncbi.nlm.nih.gov/Blast.cgi?CMD=Get&RID=Z423E71C015
https://blast.ncbi.nlm.nih.gov/Blast.cgi?CMD=Get&RID=Z42AKVWC015
The first job tells that, Sequence A belongs to botulinum neurotoxin type F [Clostridium baratii] which means that similar domain with the above botulinum neurotoxin>
Second job tells me that SequenceB belongs to unknown protein [Clostridium botulinum NCTC 2916]. So that means it is not in the same domain with the above botulinum neurotoxin right?

Looks like you've got two identical links for your blast searches

Ups. I edited the question.

What's a 'domain' here? Are you talking about protein domains? Domains of life? As in animals, plants, bacteria etc? What proteins were you given? You show DNA sequences here, not proteins. Please [edit] your question and clarify.

A protein domain is a conserved part of a given protein sequence and (tertiary) structure that can evolve, function, and exist independently of the rest of the protein chain.
This is different from global protein similarity, so similarity to neurotoxin type F doesn't necessarily mean that sequence A has similar domain to your known toxin gene. Your blast results do include conserved domain information though.

The sequence B has no conserved domain, but sequenceA and first sequence (toxin gene) has one. When I click to the conserved domain information, they show the same result. That means they are in the similar domain, right?

actually, when I click "search for similar protein architectures" in Sequence A, I see the botulinum neurotoxin type E which is the protein in the original sequence.

Sounds like you have your answer


text:
IGV jump to rsID location
Is it possible to load a list of rsID identifiers (e.g. rs66823261) to IGV to then be able to jump to the location of each of them when copy+pasting the rsID on the locus text box?
For human and other genomes, there is currently the feature of typing in a gene name, then jumping to the gene location.
I wonder if there is a way of achieving the same for rsIDs.
EDIT: thanks for the answers so far. I downloaded the Infinium GSA bed file from Illumina's website (https://support.illumina.com/downloads/infinium-global-screening-array-v1-0-support-files.html), loaded on an IGV session, which asked to create an index. I did a liftover of the file to hg38 and sorted with bedtools sort. Then I went to the text box, introduced one of the rs ids:
rs6569648
Got a positive result. Thank you all!

Yes, if you load the list you are interested in as a non-indexed "BED" file then you can jump to an rsID or any other name in the BED file. You'll need a targeted list, dbSNP for example is much too large to load in this way.

Great, thank. Is there an hg38 bed file with the rsIDs of one of the 1M SNP arrays out in the market?

you can use mysql ucsc to get the position of the SNP and then ask IGV to take some picures: see https://software.broadinstitute.org/software/igv/batch
$ echo -e 'rs66823261\nrs25' | awk '{printf("select concat(chrom,\":\",chromStart-10,\"-\",chromEnd+10) from  snp150 where name=\"%s\";\n",$1);}' | mysql --user=genome --host=genome-mysql.soe.ucsc.edu -A -D hg38 -N | awk '{printf("goto %s\nsnapshot\n",$1);}'

goto chr8:170681-170702
snapshot
goto chr7:11584131-11584152
snapshot

If you want all rsIDs (complete dbSnp database) then Pierre's answer is the best option currently.    Any bed file with all the rsIDs would be too big to load.  For a targeted list of, say, 100,000 or fewer snps you can use a bed file.
BTW igv.js already does the equivalent of Pierre's suggestion, this will be ported to IGV desktop in the near future.


text:
error in heatmap using R
I am learning about generating heatmaps using R. For this example I am using randomly generated data to test the process, but I encountered what appears to be an error or perhaps there is something am missing. Specifically, the values of Npy2r in groups A & B are the same, yet the heat map colors them differently. I am trying to understand why. Am including the R code I am using to generate the heat map
heatmap.2(x, Rowv=F, Colv = F, scale="col", col=greenred(3), key=TRUE, symkey=FALSE,
      density.info="none",cexRow=1,cexCol=1,margins=c(6,11),  trace="none",srtCol=45)

You have scale="col" in the code. What you are plotting is the z-score calculated based on the value distribution per column. Try changing it to "row" or "none".


text:
Mutations in sequences relative to reference (histogram)
I have a reference sequence of amino acids that I compare other sequences to using Bio.pairwise2.align.globalmx of Biopython package (one to many compare - one by one [is there a better way?])
It's works fine but only returns strings of the alignment with very little metadata.
Now I want to extract the indices of the amino acids that are different than the reference sequence and plot it in some kind of histogram according to protein features:
For example:

pos 1-36: 3 mutations.
pos 37-202: 4 mutations.
pos 203: 1 mutation
pos 204-660: 0 mutations

Is there any builtin package to do so?
How can I get the data of diff between 2 sequences in python? I mean which AA changed to which AA and at each position. Bio.pairwise2.align.globalmx returns some data structures of strings (see below) that can be printed but does not contain the locations in a convenient way.
LLGHKVNPFKALYNRIDMK-QVEALVEASKEE
||||||||||||||||||||||||||||||||
LLGHKVNPFKALYNRIDM-RQVEALVEASKEE

Can you post a minimal, working, and complete example of what you have gotten to work so far?

Lets start with more basic question:
How can I get the data of diff between 2 sequences in python? I mean which AA changed to which AA and at each position.
Bio.pairwise2.align.globalmx returns some data structures of strings (see below) that can be printed but does not contain the locations in a convenient way.
LLGHKVNPFKALYNRIDMK-QVEALVEASKEE
||||||||||||||||||||||||||||||||
LLGHKVNPFKALYNRIDM-RQVEALVEASKEE

The comments section is not a good place to discuss other questions. Probably best to open a new post.

Do you just want a graphical representation of the diverse regions of a multiple sequence alignment?

I'm not sure if multiple sequence alignment is what I need because I have one reference sequence that I compare all the others to. I need one-to-many comparison (which now I do by many one-to-one comparisons). Correct me if I'm wrong.

@JoeHealey, what if I wanted to do a graphical representation of the diverse regions of a multiple sequence alignment in python. Is there a package for that?

You can use AlignBuddy to process alignments in GenBank format, which allows you to keep the metadata. From the command line:
$: alignbuddy input.gb --generate_alignment mafft | \
   alignbuddy --extract_feature_sequences 'PROT_FEAT_1' | \
   alignbuddy --percent_id

If you want to do a large one-to-many comparison and loop through a bunch of features, I'd recommend importing AlignBuddy into a custom python script and calling the relevant functions directly. You can then collect the relevant data and build a histogram with matplotlib.

In answer to the comments, you can do multiple pair wise alignments and then concatenate them to visualise, if MSA is not the way to go. I have a script you can find here which will draw a graph of the Shannon entropy (uncertainty) in a column. If you create a gapped FASTA alignment for each of your pairs, and concatenate it in to a multifasta, the script should still work with it. Though it may be noisy - best just to try it.
JalView might be an option for a one-stop GUI tool, but I’m not sure if it’ll require an MSA instead of multiple pair wise.


text:
Does picard markduplicate toggle PCR duplicate samflag
I have a RNA-seq bam file and there are few reads that are puzzling me.
According to the bam header, this bam file is sorted by coordinate, created using tophat and markduplicate step is not done. But some reads are marked for being duplicate in the samflag. What is worse is when I run picard markduplicate, these reads' pcr duplicate flag is toggled, marking them not a duplicate. Also I manually found the duplicate of this read (identical reads with same start positions and mates start position) so initial marking looks true.
So my questions are:
Any idea why would this happen?
Does Tophat mark reads for being duplicate? ( I don't think so)
And does picard markduplicate toggle the duplicate flag if the reads are already marked as being duplicates?
Here is how the read looks before and after mark duplicate step.
Before:
C0RTF 1187 17 7579880 255 61M10754N40M = 7579927 10902 CTC...
0UNP1 163 17 7579880 255 61M10754N40M = 7579927 10902 CTC...

After Markduplicate:
C0RTF 163 17 7579880 255 61M10754N40M = 7579927 10902 CTC...
0UNP1 163 17 7579880 255 61M10754N40M = 7579927 10902 CTC...

And does picard markduplicate toggle the duplicate flag if the reads are already marked as being duplicates? yes: https://github.com/broadinstitute/picard/blob/master/src/main/java/picard/sam/markduplicates/MarkDuplicates.java#L342

@Pierre I know no Java, but that looks like it will unset the duplicate flag if isDuplicate is false. Not if it is already true. Am I missing something?

@terdon: if it's duplicate the flag is set https://github.com/broadinstitute/picard/blob/master/src/main/java/picard/sam/markduplicates/MarkDuplicates.java#L330 , else it is unset https://github.com/broadinstitute/picard/blob/master/src/main/java/picard/sam/markduplicates/MarkDuplicates.java#L342

@Pierre yes, exactly. If I understand correctly, the OP has a file with the duplicate flags already set but, when passed through markDuplicates, the flags that were set are now toggled: they're unset, so they are no longer marked as duplicates. The code you show would suggest that if the flag is already set, nothing will be changed. Unless picard doesn't read the flags and only does its own calculations to determine the dupes and, in this case, disagreed with whatever produced the original markings.

Yes, if MarkDuplicates encounters a pair that's marked as a duplicate that it considers (for whatever reason) to not be a duplicate then it will unset the duplicate mark. You can test this yourself by making a small BAM file either with or without duplicate entries but some marked regardless. In the case of actual duplicates, the duplicate flag will be maintained. In the case where there are no actual duplicates then the flag will be removed. At least that's the case in version 2.14.


text:
How to access organism IDS, taxon nodes from UniProt?
I try to find organism IDs used by UniProt (e.g. 9606 for "Homo sapiens") and information about the 'taxonomic kingdom' this organism belongs to. The only source on UniProt I could find is this file which contains all required information. So in the Python world it would be a dictionary of the form:
{org_id: {'name': 'org_name', 'kingdom': 'org_kingdom'}}

or a dataframe.
However, the file from above is not particular nice to parse due to its inconsistent column format (not impossible to parse though), so I am wondering whether there is an alternative source to retrieve this information or whether there is maybe even a webservice available that takes care of this parsing.

That file is actually trivial to parse. You only care about lines that have N=$species, so you can simply do:
sed -En 's/.* ([0-9]+): N=(.*)/\1\t\2/p' speclist

That will return a tab separated list of taxID and species:
$ sed -En 's/.* ([0-9]+): N=(.*)/\1\t\2/p' speclist | head
648330  Aedes albopictus densovirus (isolate Boublik/1994)
10804   Adeno-associated virus 2
648242  Adeno-associated virus 2 (isolate Srivastava/1982)
118452  Abacion magnum
72259   Abaeis nicippe
102642  Abax parallelepipedus
392897  Abalistes stellaris
75332   Abbottina rivularis
515833  Abdopus aculeatus
56673   Antarctic bacterium DS2-3R

So to get human, for example, you could do:
$ sed -En 's/.* ([0-9]+): N=(.*)/\1\t\2/p' speclist | grep sapiens
63221   Homo sapiens neanderthalensis
9606    Homo sapiens

And, for exact matches:
$ sed -En 's/.* ([0-9]+): N=(.*)/\1\t\2/p' speclist |
    awk -F'\t' '$2=="Homo sapiens"'
9606    Homo sapiens

If you also want to keep the taxon node (kingdom), you can use:
sed -En 's/.* ([A-Z])  *([0-9]+): N=(.*)/\1\t\2\t\3/p' speclist

Which prints:
$ sed -En 's/.* ([A-Z])  *([0-9]+): N=(.*)/\1\t\2\t\3/p' speclist  | head
V   648330  Aedes albopictus densovirus (isolate Boublik/1994)
V   10804   Adeno-associated virus 2
V   648242  Adeno-associated virus 2 (isolate Srivastava/1982)
E   118452  Abacion magnum
E   72259   Abaeis nicippe
E   102642  Abax parallelepipedus
E   392897  Abalistes stellaris
E   75332   Abbottina rivularis
E   515833  Abdopus aculeatus
B   56673   Antarctic bacterium DS2-3R

The sed command might appear a little daunting if you don't know sed, but it's actually quite simple. The -E flag enables extended regular expressions which give us a simpler syntax. The -n suppresses normal output, so that sed won't actually print anything unless explicitly told to.
The s/from/to/ is the substitution operator and will replace from with to. Here, we are replacing zero or more characters (.*) until a space followed by one or more numbers ([0-9]+), a :, another space, then the string N= and then everything else (.*). The parentheses let us capture the matched strings, so that we can then refer to the 1st captured string as \1, the second as \2 etc. We are therefore replacing the entire line (because the regex matches the entire line) with the 1st pattern captured (the taxID, \1), a tab (\t) and the species name (\2). Finally, the p at the end means "print the resulting line if the substitution was successful", so only lines matching the desired pattern are actually printed.

Thanks, that helps already. Is there a quick way to also extract the taxonomic kingdom and store it in the third column?

@Cleb sure, see updated answer. Next time, please explain all these details in your question (I mean tell us exactly what information you want to extract) so we can give you a working answer directly.

These are NCBI taxonomy ID's rather than anything specific to UniProt. One way to parse them would be using BioPython via Bio.Entrez, for example:
from Bio import Entrez

def get_taxon_name_from_NCBI_id(id_):
    Entrez.email = 'name@provider.com'
    handle = Entrez.efetch(db='taxonomy', id= id_, retmode='xml')
    records = Entrez.read(handle)[0]
    return records['ScientificName'], records['TaxId']

print(get_taxon_name_from_NCBI_id(9606))
# ('Homo sapiens', '9606')

Of course you could query Entrez directly:
https://eutils.ncbi.nlm.nih.gov/entrez/eutils/efetch.fcgi?db=taxonomy&id=9606
Or even more basically do a manual search of NCBI taxonomy:
https://www.ncbi.nlm.nih.gov/taxonomy/?term=9606

I was not precise enough in my question: What I am after are the IDs, so basically what you use as input and there associated name and taxonomic kingdom. Ideally it would be a dictionary of the form {org_id: {'name': 'org_name', 'kingdom': 'org_kingdom'}} but I guess I can easily do this directly from NCBI; thanks for making me aware of this!