Hot answers tagged

147

Good observation! The 3' poly(A) tail is actually a very common feature of positive-strand RNA viruses, including coronaviruses and picornaviruses. For coronaviruses in particular, we know that the poly(A) tail is required for replication, functioning in conjunction with the 3' untranslated region (UTR) as a cis-acting signal for negative strand synthesis ...


51

Let’s start with what they have in common: All three formats store sequence data, and sequence metadata. Furthermore, all three formats are text-based. However, beyond that all three formats are different and serve different purposes. Let’s start with the simplest format: FASTA FASTA stores a variable number of sequence records, and for each record it ...


44

That is the correct sequence for 2019-nCov. Coronavirus is of course an RNA virus and in fact, to my knowledge, every RNA virus in Genbank is present as cDNA (AGCT, i.e. thydmine) and not RNA (AGCU, i.e. uracil). The reason is simple, we never sequence directly from RNA because RNA is too unstable and easily degraded by RNase. Instead the genome is reverse ...


31

This question is quite general, so I'm going to attempt to tie it back to bioinformatics. Background The tree for the current coronavirus is here, showing it is closely related to bat-coronavirus and in particular SARS. Question The bioinformatics question for the current coronavirus is why this virus appears to be able to infect humans and transmit to ...


26

What does this soft masking actually mean? A lot of the sequence in genomes are repetitive. Human genome, for example, has (at least) two-third repetitive elements.[1]. These repetitive elements are soft-masked by converting the upper case letters to lower case. An important use-case of these soft-masked bases will be in homology searches: An atatatatatat ...


24

For FASTQ: seqtk fqchk in.fq | head -2 It gives you percentage of "N" bases, not the exact count, though. For FASTA: seqtk comp in.fa | awk '{x+=$9}END{print x}' This command line also works with FASTQ, but it will be slower as awk is slow. EDIT: ok, based on @BaCH's reminder, here we go (you need kseq.h to compile): // to compile: gcc -O2 -o count-N ...


24

The scenarios are impossible and would be laughable if they were not so serious. The evidence is in the phylogenetic trees. Its a bit like a crime scene when the forensics team investigate. We've done enough crime-scenes often going to the site, collecting the pathogen, sequencing and then analysis - (usually neglected diseases) without any associated ...


20

Some of the other answers here seem quite good; at the same time I think the core answer to the OP's question is maybe a bit hard to tease out of them, so I'd like to try to state it more plainly. It's worth noting that a truly complete answer to this question seems to be beyond current research, but any kind of "Why?" is inevitably a hard or even impossible ...


19

If you want something quick and dirty you could rapidly index the FASTA with samtools faidx and then put the lengths column through R (other languages are available) on the command line. samtools faidx $fasta cut -f2 $fasta.fai | Rscript -e 'data <- as.numeric (readLines ("stdin")); summary(data); hist(data)' This outputs a statistical summary, and ...


18

5 hours and no benchmarks posted? I am sorely disappointed. I'll restrict the comparison to just be fasta files, since fastq will end up being the same. So far, the contenders are: R with the ShortRead package (even if not the fastest, certainly a super convenient method). A pipeline of grep -v "^>" | tr -cd A | wc -c A pipeline of grep -v "^>" | ...


17

In a nutshell, FASTA file format is a DNA sequence format for specifying or representing DNA sequences and was first described by Pearson (Pearson,W.R. and Lipman,D.J. (1988) Improved tools for biological sequence comparison. Proc. Natl Acad. Sci. USA, 85, 2444–2448) FASTQ is another DNA sequence file format that extends the FASTA format with the ability ...


16

Most sequencing experiments, be it Illumina-based next-generation-sequencing or Sanger sequencing uses DNA as template, not RNA. Even if this virus is RNA-based it would be reverse-transcribed prior to any sequencing experiment. Therefore the output is DNA and this is what NCBI provides here.


14

You can do this easily with bioawk, which is a version of awk with added features facilitating bioinformatics: bioawk -c fastx '{print $name"\t0\t"length($seq)}' test.fa -c fastx tells the program that the data should be parsed as fasta or fastq format. This makes the $name and $seq variables available in the awk commands.


13

Not an expert, but some searching on eukaryotic positive-strand RNA viruses seems to show that polyadenylation is not uncommon. For example, Steil, et al., 2010.


12

Statistics for nanopore reads are tricky because of the huge range of read lengths that can be present in a single run. I have found that the best way to display lengths is by using a log scale on both the x axis (length) and the y axis (sequenced bases, or counts, depending on preference). I have written my own scripts for doing this: one for generating ...


12

There's rarely a good reason to use a hard-masked genome (sometimes for blast, but that's it). For that reason, we use soft-masked genomes, which only have the benefit of showing roughly where repeats are (we never make use of this for our *-seq experiments, but it's there in case we ever want to). For primary vs. toplevel, very few aligners can properly ...


11

It's good practice to have your FASTA indexed, so you can leverage the .fai you are likely to already have. If not, you can just generate the index with samtools and use some awk to make your BED: samtools faidx $fasta awk 'BEGIN {FS="\t"}; {print $1 FS "0" FS $2}' $fasta.fai > $fasta.bed This will maintain tab separation but you can drop the BEGIN ...


11

The use of lower/upper case letters and N/n letters in genomes sequences is not completely standardised and you should always check the specification of the resource you are using. Lower case letters are most commonly used to represent “soft-masked sequences”, a convention popularised by RepeatMasker, where interspersed repeats (which covers transposons, ...


11

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,...


10

If I use SeqIO.parse(filehandle, 'fasta') to parse a FASTA file, then it will return a SeqRecord object where the id and name are the first word (everything before the first whitespace) of the line beginning with > and the description is the complete line (all not including the initial >). (This behaviour can overruled by providing a custom title2ids ...


10

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))


9

I think that this should be pretty fast: FASTA: grep -v "^>" seqs.fa | tr -cd N | wc -c FASTQ: sed -n '1d;N;N;N;P;d' seqs.fq | tr -cd N | wc -c See this answer on SO about how to count characters in BASH using different approaches.


9

While there's nothing stopping anyone from doing that with the FASTA format (after all, it's just a text file with '>' defining header lines), I don't know of any software that would support such a file structure. At best, it would interpret the nucleotide sequences as protein sequences (A/C/G/T are all valid 1-letter protein codes). A better question to ...


9

Why not use sed instead? sed -e 's/chr_I/I/' -e 's/chr_V/V/' -e 's/chr_X/X/' mySequence.fasta > mySeq.fasta Or even simpler: sed 's/chr_//' mySequence.fasta > mySeq.fasta


9

If you want to avoid using extra libraries for any reason, you can just use a simple Python script (version 3.6 and above) to do this: fr = open("dup_test.fasta", "r") fw = open("dup_edited.fasta", "w") seq_dict = {} curr_header = '' for line in fr: line = line.strip() if line[0] == '>': if line not ...


8

FASTQ As it was pointed out, fastq can be complicated. But in a simple case when you have four lines per record, one possible solution in bash is: sed -n '2~4p' seqs.fastq | grep -io N | wc -l sed -n '2~4p' will print every fourth line grep -o N will output a line with N for every matching symbol wc -l will count the lines I suspect this python approach ...


8

Generally, you should use the soft-masked or unmasked primary assembly. Cross-species whole-genome aligners, especially older ones, do need to know soft-masked regions; otherwise they can be impractically slow for mammalian genomes. Modern read aligners are designed to work with repeats efficiently and therefore they don't need to see the soft mask. For ...


8

It's not a fasta file, but: > m sample1 sample2 sample3 sample4 fliI 1 1 1 1 patB_1 1 1 1 1 pgpA 1 1 1 1 osmB 1 1 1 1 cspA 1 0 1 1 > # Collapse to labeled strings > blah = apply(m, 2, function(x) paste(x, collapse=''...


8

The description field in the SeqRecord object has the information you are looking for: >> from Bio import SeqIO >> s = SeqIO.read('genome.fasta', 'fasta') # single sequence fasta file >> s.id k119_5 >> s.description k119_5 flag=0 multi=141.0706 len=473 Edit: as an aside, if you write a SeqRecord object using SeqIO's write method, ...


8

This looks like a bug in makeblastdb. Removing the final | | | from your sequence's description makes it work: >K02545|TRBJ1-1*01|Homo sapiens|F|J-REGION|749..796|48 nt|3| | | |15 AA|15+0=15 So does removing everything after the 1st space: >K02545|TRBJ1-1*01|Homo sapiens|F|J-REGION|749..796|48 There is nothing in the FASTA definition that would ...


Only top voted, non community-wiki answers of a minimum length are eligible