Tag Info

22

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

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

11

You can use the Jellyfish software to calculate the k-mer profiles up to length 31. From the instructions in the user guide: The basic command to count all k-mers is as follows: jellyfish count -m 21 -s 100M -t 10 -C reads.fasta To compute the histogram of the k-mer occurrences, use the histo subcommand (see section 3.1): jellyfish histo mer_counts.jf ...

10

You can also use R. I give you an example of only chr1 and only kmer=4. library(BSgenome.Hsapiens.UCSC.hg38) library(Biostrings) genome <- BSgenome.Hsapiens.UCSC.hg38 kmers <- oligonucleotideFrequency(genome$chr1, 4) kmers m <- as.matrix(kmers) m[order(m),] 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

For a quick estimate you’re making it more complicated than necessary. The theoretical average coverage is $\frac{n \cdot \hat l}{N}$ where $n$ is the number of reads, $\hat l$ is the average read length and $N$ is the genome size. samtools idxstats gives you the chromosome lengths and number of mapped reads in one convenient list. Putting this together, ...

7

I would also recommend to take a look at pomegranate, a nice Python package for probabilistic graphical models. It includes solvers for HMMs and much more. Under the hood it uses Cythonised code, so it's also quite fast.

6

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

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

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

5

What you are looking for is SNP annotation. If you have the chromosome:position reference and alternate alleles for your SNPs of interest, it can be as simple as uploading them to the variant effect predictor. http://grch37.ensembl.org/Homo_sapiens/Tools/VEP This will give you the predicted protein change and novelty of the variant with respect to known ...

5

awk awk '{L=length($1);for(i=1;i<=L;i++) {B=substr($1,i,1);T[i][B]++;}} END{for(BI=0;BI<4;BI++) {B=(BI==0?"A":(BI==1?"C":(BI==2?"G":"T")));printf("%s",B); for(i in T) {tot=0.0;for(B2 in T[i]){tot+=T[i][B2];}printf("\t%0.2f",(T[i][B]/tot));} printf("\n");}}' input.txt A 0.13 0.40 0.13 0.00 0.00 0.40 0.47 0.07 0.40 C 0.40 ...

5

There are certainly software libraries for working with HMMs. For a general-purpose implementation in C++, take a look at the SeqAn HMM algorithms. For your purposes, i.e. “computing … the most likely hidden sequence given the sequence of observed states”, you’d invoke viterbiAlgorithm with your observed sequence and the HMM graph. More fundamentally I ...

5

Now there's another option, the ensembl_rest module, a thin wrapper around the Ensembl REST API to simplify its usage and make it more pythonic. You can find the documentation here. To clarify things, I'm the creator and maintainer, but still think it's a legitimate alternative.

5

Panel (a) shows a sequence gap-free alignment. Each row corresponds to a contiguous 15 base pair sequence of DNA (e.g. row 1 could be a human sequence, row 2 could be the equivalent mouse sequence etc.). Each column corresponds to a particular position/residue (e.g. column 1 is the 1st position of the alignment etc.). Panel (b) shows the sequence logo. As ...

5

I can only speak of drug design (and even then I am terrible at turning down the jargon). In the case of drug design, this is pretty much plan C. Namely, none of compounds that entered clinical trial at the start of the year work (let's call this plan A) and none of the vaccines that are entering now clinical trial work (let's call this plan B although ...

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

If I remember correctly Ewan Birney's Dynamite (a compiler-compiler) as presented at ISMB 1997 had this functionality, there is also some code here on GitHub https://github.com/birney/wise3 which at least mentions Dynamite. Suspect Ewan is too busy these days to work on this, although he has tweeted about blowing dust of his old Dynamite, sorry Dynamite ...

4

If I understood correctly your question, you want to mask those regions in a (FASTA?) genome. I think you could identify those regions using mummer and mask them using bedtools. # align genome against itself nucmer --maxmatch --nosimplify genome.fasta genome.fasta # select repeats and convert the corrdinates to bed format show-coords -r -T -H out.delta | ...

4

Here is one example of how to do this with a bit of python. Alternatively one could create strings of each column and using letterFrequency() from the Biostrings package. #Make a list of hashes hl = [] for i in range(9): hl.append({'A': 0, 'C': 0, 'G': 0, 'T': 0}) f = open("foo.txt") # CHANGE ME nLines = 0 for line in f: for idx, c in enumerate(...

4

Grabbing chromosomes for hg38: $wget ftp://hgdownload.cse.ucsc.edu/goldenPath/hg38/chromosomes/*.fa.gz$ for fn in ls *.fa.gz; do gunzip fn; done Via kmer-counter and Python, here's how to search for kmers of length 7 from chromosome chrY: #!/usr/bin/env python ... 4 You can drastically simplify your effort because FASTQ files are simple text files; therefore, standard UNIX test file tools work: cat refs/*.fastq > combined.fastq And this even works with .fastq.gz files — no unzipping required. Apart from this, your R code contains several errors. For instance, your variable names change throughout, and you’re using ... 4 Here are a couple books I'd recommend: Dan Gusfield's Algorithms on Strings, Trees, and Sequences is a deep and wide treatment of aligning, searching, and processing strings, trees, and sequences. Warren Ewens' Statistical Methods in Bioinformatics devotes a chapter to BLAST and the math underneath it. Edit - Another book that may be useful for some ... 4 cellranger mkfastq is not necessary anymore. It used to be that the cellranger software wanted the reads to be interleaved, and you could use cellranger to do that for you if you couldn't do it yourself. Newer versions of cellranger will take the fastq files just like Illumina's bcl2fastq makes them. Cellranger count aligns the reads, filters away ... 4 The correct API for Ensembl is the Ensembl REST API which is updated and maintained by Ensembl, and language agnostic. 4 The RAM requirements for some bioinformatics analyses like assembly can be quite high (in the hundreds of Gigabytes). My recommendation is to get a fast laptop. Something with an i7 quad-core processor, 16 GB of RAM, and 1TB of storage should do. Do your analyses on this laptop. If you run out of RAM or storage, then use a cloud service to get access to ... 4 While I haven't found a way to limit the results to the canonical transcript only, you can get a list of genes, transcripts and their CDS lengths using Ensemble's BioMart. I have already set it up for you, you can see the results, and modify them, here (click on the "Results" link if you don't see them). Essentially, you just need to go to BioMart, and ... 4 You can use al$-order Markov chain. Here is the procedure: Count$l$-mers in your genome. For small genomes, you can do that in Python. For large genomes, you may need jellyfish or KMC3. Draw a$l$-mer randomly based on the distribution of$l$-mers. This generates the first$l$bases. Let$s$be the last$(l-1)\$ subsequence from the generated sequence. ...

4

Why do you think the entropy of 0 is incorrect? It intuitively makes sense, as there is no uncertainty about the base at position 3, and thus there is no entropy. However, what is plotted in a sequence logo isn't the entropy, but rather a measure of the "decrease in uncertainty" as the sequence is aligned. This is calculated by taking the entropy at this ...

4

Yes, if the insert size is smaller than the read sizes, this would happen. For some applications, for example SNP detection for molecular diagnostics purposes, this approach is used.

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