Hot answers tagged

30

The earliest mention of the 30x paradigm I could find is in the original Illumina whole-genome sequencing paper: Bentley, 2008. Specifically, in Figure 5, they show that most SNPs have been found, and that there are few uncovered/uncalled bases by the time you reach 30x: These days, 30x is still a common standard, but large-scale germline sequencing ...


23

I don't know if it's the fastest, but the following provides an approximately 10x speed up over your functions: import string tab = string.maketrans("ACTG", "TGAC") def reverse_complement_table(seq): return seq.translate(tab)[::-1] The thing with hashing is that it adds a good bit of overhead for a replacement set this small. For what it's worth, I ...


21

Thanks to Manu Tamminen for this solution: echo ACCTTGAAA | tr ACGTacgt TGCAtgca | rev


21

All humans have some differences in their DNA, but there's far more that is shared. On average the difference between humans is only about one thousandth of their full DNA, which means we're about 99.9% the same. These differences aren't distributed fully randomly, but are often because of specific gene alternatives. (Random mutations do occur, but they are ...


18

The HGP developed the first "reference" human genome - a genome that other genomes could be compared to, and was actually a composite of multiple human genome sequences. The standard human reference genome is actually continually updated with major and minor revisions, a bit like software. The latest major version is called GRCh38, was released in 2013, ...


16

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


9

Here's a Cython approach that might suggest a generic approach to speeding up Python work. If you're manipulating (ASCII) character strings and performance is a design consideration, then C or Perl are probably preferred options to Python. In any case, this Cython test uses Python 3.6.3: $ python --version Python 3.6.3 :: Anaconda custom (64-bit) To ...


8

While the quality of the reference human assembly keeps improving, there are still misassemblies in it. A common problem is recent segmental duplications are occasionally collapsed into one sequence in the reference. Another issue is that the centromeric sequences in the reference are computationally generated, which are probably different from real ...


8

using a http request. if there is a DAS server, you can always use this protocol to download the xml -> fasta. see https://www.biostars.org/p/56/ $ curl -s "http://genome.ucsc.edu/cgi-bin/das/hg19/dna?segment=chrM:100,200" | xmllint --xpath '/DASDNA/SEQUENCE/DNA/text()' - | tr -d '\n' ...


8

DNA is only the substance that makes up the genome. DNA can come with many forms. For example, the fragments you get by PCR in the lab are DNA. A bacterium also has DNA molecules in it. You may also find DNA fragments in many other places, like in our blood, or at a crime scene. Genome, however, is a very specific term. It means all the heritable ...


8

TLDR; For the most modern CHM13 assembly (v1.1; Complete T2T reconstruction of a human genome; https://github.com/marbl/CHM13), the best compression without any use of a reference would be between 624 MB and 641 MB. Upper estimates of genome compressibility can be made using standard compressors (see, e.g., this comparison). In my experience, xz can achieve ...


7

The most reliable and simplest way is probably using Biopython: from Bio.Seq import Seq def my_reverse_complement(seq): return Seq(seq).reverse_complement() print(my_reverse_complement('ATGCGTA')) # TACGCAT As Devon has already said here using Biopython isn't as fast as the naive Python solution, and I also tested that shown here with ipython. ...


7

Solexa Inc. sequenced NA12878 chrX to ~30x in early 2007, which later became part of Bentley (2008). This, I believe, was the first time that 30x showed up. I don't recall they had a particular reason for that. Figure 5 in the published paper was just aftermath. It does not really explain why not 25x or 35x, given that the curves between 25x and 35x in that ...


7

There are a few quick'n'dirty ways depending on the type of data. In any case you want to align your files to a reference genome and then check the distribution of reads, either on a genome browser or with tools such as RSEQC which calcualtes the fraction of reads aligning to exon, intron, intergenic etc. RNA-seq, if you use a standard aligner such as ...


6

30 times coverage is not unique to this problem, but number 30 has its empirical role in statistics: In statistical analysis, the rule of three states that if a certain event did not occur in a sample with n subjects, the interval from 0 to 3/n is a 95% confidence interval for the rate of occurrences in the population. When n is greater than 30, this is a ...


6

Not as far as I am aware. The Ray assembler used to (and possibly still does) store the kmers as FASTA files where the header was the count of the sequence, which I thought was a pretty neat bastardisation of the FASTA file format. It looks like this format is also used by Jellyfish when reporting kmer frequencies by the dump command (but its default output ...


6

I was quite curious when I saw that as well. I've spent more time than I'd care to admit trying the sort of things you have. I've gotten it decoded now, but I can't really claim any sort of victory, as the problem was pre-solved. After struggling through some of the same experiments you did, I decided to take a closer look at their explainer video here: ...


5

Another python extension but without cython. The source code is available at the bottom of this answer or from this gist. On Mac with Python3: string 0.38166s total, 1310067.4 strings per second 84.4% increase over baseline seqpy 0.20524s total, 2436182.7 strings per second 91.6% increase over baseline On Linux with Python2 (seqpy is the ...


5

Based on my experience reviewing k-mer counting software and as a core contributor to the khmer project, I can confidently say that there is no widely used standard format.


5

Early MinION sequencing runs had forward and reverse DNA templates joined together by a hairpin adapter, so that the sequencer would read both strands from the same template. The consensus sequence that was generated by combining the base calls of opposite strands is referred to as a 2D read. Due to a dispute with Pacific Biosciences, and for other reasons ...


5

It might be used differently in different contexts, but generally speaking, in my world - allelic imbalance is when there's a difference in the level of gene expression from different alleles, usually through genetic (e.g. a variant in a promoter) or epigenetic mechanisms (e.g. one copy silenced, as in imprinted regions). There are a few blog posts here ...


5

Reverse complement FASTA/Q: seqtk seq -r in.fa > out.fa https://github.com/lh3/seqtk


5

If I'm understanding you correctly, by "types of DNA methylation" you mean "nucleotide contexts where DNA methylation occurs". This is going to be a function of the methyltransferase proteins involved in the methylation process, and this largely depends on the organism, and even cell type being studied. Generally, what has been observed is: The phenomenon ...


5

Although there is not a unique nucleotide sequence that translates to a given protein, one can list all the possible DNA sequences that do translate to that protein. An online tool that does just that is Backtranambig, from EMBOSS. It produces a DNA sequence representing all the nucleotide sequences matching the input protein, using IUPAC ambiguity codes.


5

There's a benchmark that compared dozens of compressors on DNA data, including human genome: http://kirr.dyndns.org/sequence-compression-benchmark/ In particular, this table shows file size of compressed human genome, produced by the strongest setting of each compressor: http://kirr.dyndns.org/sequence-compression-benchmark/?d=Homo+sapiens+GCA_000001405.28+%...


5

To get around limitations in using Venn diagrams for set overlaps, when there are more than three sets, back around 2013 or so I created something I called an Eulergrid plot (example at the bottom of the page), which an UpsetR plot appears to recapitulate, today. The Eulergrid code I wrote was in a mix of Perl and R; the UpsetR plot code uses R. There ...


5

I would just blast it. When blasting locally, you need to first make a database from your genome, so assuming you got the command-line version of blast installed you can do something like makeblastdb -in my_study_genome.fa -dbtype nucl blastn -max_target_seqs 10 -db my_study_genome.fa -evalue 1e-10 -outfmt 6 -query my_downloaded_gene_of_interest.fasta -out ...


4

Depending on the coverage of your data and the complexity of the genome, you could either reassemble the genome de novo or run a reference guided (or reference assisted) assembly. It sounds like you're leaning more towards the latter. There are a couple of reference-guided assembly tools available: AlignGraph and Ragout. These may or may not be appropriate ...


4

One approach to this is to use whatever data you have to iteratively update the reference genome. You can keep chain files along the way so you can convert coordinates (e.g. in gff files) from the original reference to your new pseudoreference. A simple approach might be: Align new data to existing reference Call variants (e.g. samtools mpileup, GATK, or ...


4

There are more R packages available that calculate GC content, for example Ape's GC.content() function. For example: > library(ape) > data(woodmouse) > GC.content(woodmouse) [1] 0.3873347 With a sliding window is in Biostrings package. > library(Biostrings) > DNA <- DNAString("ACTGAAACCGTGGCAGTTTGAC") > letterFrequencyInSlidingView(...


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