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


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


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


16

See IUPAC codes: So, as you can see above, K means "Either G or T".


14

First of all, if you want to understand mapping quality (mapQ), ignore RNA-seq mappers. They often produce misleading mapQ because mapQ is not important to RNA-seq anyway. Strictly speaking, you have two questions, one in the title: the meaning of mapQ; and the other in a comment: how mapQ is computed. On the meaning, mapQ is nearly the same as baseQ – the ...


14

Bowtie2 is no longer the fastest aligner. Salmon and Kallisto are much faster, but have been designed to optimise RNASeq mapping. Their speed is gained from avoiding a strict base-to-base alignment, but they can output mostly-aligned reads (i.e. position-only, without local alignment) as pseudo-alignments. See here for more details. Both Kallisto and Salmon ...


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

The seed is the subset of a read used in the first step of an alignment. Many aligners work by a seed-and-extend model, wherein they first find all regions matching the "seed" and then extend the alignment around that allowing mistmatches and indels until it either gives up (and therefore uses a different seed) or finds a sufficiently good alignment.


11

There are three possible problems that come to mind. Masking Blast will mask low complexity regions by default. Since your sequence is nothing but Gs, it is a safe bet that it is being masked, so no hits will be found for it. Score/e-value thresholds Another source of complication is that even if a match is found, that match will have very bad scores. Both ...


11

You are looking for the needle program from the EMBOSS suite. Available in bioconda. http://emboss.sourceforge.net/apps/release/6.6/emboss/apps/needle.html To read sequences from the commandline, you need so specify the format as asis. To get the output on the screen, you'll need -stdout and to use the default alignment parameters (gap penalty 10, ...


11

UPDATE: The article has now been withdrawn with the following note: This paper has been withdrawn by its authors. They intend to revise it in response to comments received from the research community on their technical approach and their interpretation of the results. If you have any questions, please contact the corresponding author. This is very ...


10

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: #!...


9

I'm trying to figure out why. This will be a longer read, Tl;dr at the end: Doing a match There is a good match of Q against LC074724 in 20170824: $ blastn -outfmt '6 qaccver saccver bitscore' -db 20170824 -query Q.fasta \ -task blastn -max_hsps 1 | grep LC074724 Q LC074724.1 2581 Note that the bitscore of the LC074724 match is 2581. But there is ...


9

Install biopython with conda and use the following script: #!/usr/bin/env python import sys from Bio import Align aligner = Align.PairwiseAligner() aligner.mode = "local" alignments = aligner.align(sys.argv[1], sys.argv[2]) print(alignments[0]) The output is then: .GGA-GGAGGGAG--AAGGAGGGAGGGAAGA-GGAGGGAG--AAGGAGGGAGGC .|-|-||||||||--|||-|||-|||||-||-||||||...


8

I am not sure what you mean by "fasta alignment file". If you mean a multi-sequence alignment (MSA) in the fasta format, you can't get that because SAM keeps pairwise alignments only and doesn't align inserted sequences. Even if you don't care about inserted sequences, a MSA in fasta is far to big to be practical. Alternatively, by "fasta alignment file", ...


7

SNPs are likely to be created and InDels are likely to be missed. Suppose you have a read, ACTGACTGACTGTAC and you align it to a reference sequence ACTGACTGACTGTTAAGAACGACTACGAC. If you aligned that, you would either get: ACTGACTGACTGTac (lower case denotes soft-clipping) ACTGACTGACTGTTAAGAACGACTACGAC or ACTGACTGACTGTAC (N.B., you've created some ...


7

The first column is added in order to be able to align the sequences. They might be gaps if the alignment ends up beginning in different positions than the first. Without the first empty row and column it wouldn't be possible to have gaps at the beginning of the alignment The dynamic programming is useful because you need first to calculate the score for ...


7

Insertions and deletions are the dominant error mode of long read sequencing, including nanopore sequencing. What you see is not unexpected. Things may have improved by now if you would download the raw fast5 data and repeat the basecalling. There is no need to gunzip the fastq.gz prior to alignment. Your commands for alignment look alright, except that (if ...


7

I would not look for a package for this, but instead build a small pipeline calling external tools with something like the following workflow: Cluster the ~100 sequences with CD-HIT-EST/PSI-CD-HIT or many other options Take all the sequences that form one individual cluster and build a multiple sequence alignment (MSA) with MAFFT/ClustalOmega or similar ...


7

Normally "inserts" used in the manuscript are "indels" in protein alignments, short for insertions and deletions. What I think has happened is a group investigating indels in HIV env noticed indels in 2019-nCov. Essentially I think the correlation is spurious - but I haven't test it, but the area of research in understanding indels is certainly valid and ...


6

For simple variants like SNPs it would not really be a problem to use the current genome assembly for other ethnic groups. But for more complex variants this could be indeed problematic, however not only for ethnic groups but also for individuals within the same population. Think of very complex regions such as HLA or KIR. In studies where they compare ...


6

split reads - These are read that have two or more alignments to the reference from unique region of the read. In this example a 150bp read sequenced from RNA could have base 1-75 aligning to the 3' end of exon2 and bases 76-150 aligning to the 5' end of exon3. This would be a split read because it have two alignments (exon2 and exon3) and those alignment ...


6

Found a solution, using D-Genies, worked great. Some examples from their website: Thanks to @user172818.


6

BLAT can only use one CPU. It is actually not the right tool for full-genome alignment. For "two versions" of the same species, MUMmer and minimap2 are orders of magnitude faster and probably give better alignment. EDIT (moving comment to answer): OP comments that the purpose of this alignment is for lifting over annotations using the UCSC PSL-based ...


6

One possibility is to run bcftools norm on the file. At the end it will print out a statistic about how many variants were realigned. I did this for dbSNPv151 (GRCh38) only for chromosome 22 like this: $ bcftools view -i 'TYPE="indel"' dbSNP_GRCh38.p7_b151.vcf.gz 22|bcftools norm -Ou -f GRCh38.fa > /dev/null Lines total/split/realigned/skipped: 901503/...


6

Bits are not frequencies. If a position only contains an A (position 3 for example) then you would need 2 questions (bits) to derive that information without a priori knowledge. Is it a G or C? if no > then it is a A or T Is it T? If no then it is an A Position 1 can be derived by 1 question only: is it a G or C? If no then it is an A or T In this case ...


6

Overview The central focus of the tree is to highlight the key biological concern of the new coronavirus, 2019-nCov. The key concern is the genetic similarities to SARS epidemic, and relates to the SARS receptor. SARS background SARS is endemic in bats (your BioRxiv tree partly shows that and this tree definitely shows it) and in the 2002 epidemic infected ...


6

You are right that a single molecule in a single position is either methylated or not methylated. However: First, assuming your organism of interest is diploid (or of higher ploidy) one of the chromosomes could be methylated, the other not. That would give you a level of 0.5 and can be found in imprinted regions (where the paternally inherited chromosome is ...


6

Clarification Just making sure - you have ~6000 datasets - one for each S.cerevisiae gene of interest - which are made up of homologues of each respective gene? And you want to filter out genes that are too similar? In that case, it would be helpful to know how you're defining a sequence as 'similar enough' for it to be filtered. Given they're all homologues,...


5

The linear gap penalty has a cost which is directly proportional to the length of the gap: the gap costs 1 unit per consecutive missing base. The simplest case is assigning each unit (= each missing base) a penalty of 1: g(k) = k Or more generally, if you want to increase the penalty per base to 2, 3 or n: g(k) = n * k with k = length of the gap, ...


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