Hot answers tagged

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


18

I'm not sure I'm doing it the best way, but here is an example where I read a compressed gzip fastq file and write the records in block gzip fastq: from Bio import SeqIO, bgzf # Used to convert the fastq stream into a file handle from io import StringIO from gzip import open as gzopen records = SeqIO.parse( # There is actually simpler (thanks @peterjc) ...


18

For Biopython 1.70, there is a new Seq.count_overlap() method, which includes optional start and end arguments: >>> from Bio.Seq import Seq >>> Seq('AAAA').count_overlap('AA') 3 >>> Seq('AAAA').count_overlap('AA', 1, 4) 2 This method is also implemented for the MutableSeq and UnknownSeq classes: >>> from Bio.Seq import ...


13

Taking a different tack from other answers, there's lots of tools for pipelines in Python. Note: there was a time when people would use "pipeline" to refer to a shell script. I'm talking about something more sophisticated that helps you decompose an analysis into parts and runs it robustly. Snakemake is my favourite. It's (nearly) pure Python and can ...


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


10

A simple solution for converting Seurat objects into AnnData, as described in this vignette: library(reticulate) seuratobject_ad <- Convert(from=seuratobject, to="anndata", filename="seuratobject.h5ad") Alternatively, one can use Loom, "a file format designed to store and work with single-cell RNA-seq data". In R, install LoomR: devtools::...


9

I've encountered this problem before, and used python re module to solve this problem. import re all = re.findall(r'(?=(AA))','AAAA') counts = len(all) You can get more details in this thread


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

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

The first place to start is the GFF3 specification. This is the official word on what is and is not allowed in a GFF3 file. For example, users can define arbitrary attribute keys, so long as they do not begin with an uppercase letter (these are reserved for "official" use). But your question doesn't seem to be about what is allowed, but what is commonly ...


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

The trick would be to swap the key in the dictionary to be the sequence itself. Also I would recommend using a different separator that "_" since that is what the current ids have so that you can easily separate the individual ids from the concatenated id. I used a pipe "|" in this example. Also I just manually wrote the FASTA output instead of using ...


8

To expand on my comment from yesterday. You could do this with the ETE Toolkit (I just copied one logo file rather than converting all 26 to png): from ete3 import Tree, TreeStyle, faces def mylayout(node): if node.is_leaf(): logo_face = faces.ImgFace(str.split(node.name, '.')[0] + ".png") # this doesn't seem to work with eps files. You could ...


8

you are starting snakemake the wrong way. snakemake snakemake means "Running the pipeline in Snakefile until you reach a rule named snakemake or a file name snakemake is created." If your snakefile is called snakemake and this is the pipeline you want to start, the correct syntax would be: snakemake -s snakemake


8

You're assuming 2-line FASTA format, but the sequence can (and typically does) span multiple lines. You can use Biopython to build a simple parser: from Bio import SeqIO for record in SeqIO.parse("file_name.fa", "fasta"): if "16S" in record.description: print(record.format("fasta"))


7

Splitting into multiple files and changing the IDs can be easily done: perl -pe 'if(/>/){/\[(.*?)\]\s*$/; $_="> $1\n"}' file.fa | awk '(/^>/){name=$2} {print >> name".fa"}' That assumes all your FASTA headers have [foo bar baz] as the last element of a line. It will create a file called foo.fa (the bacterium's name) with all sequences ...


7

You can convert any string hash function to a "canonical" DNA string hash function. Given a DNA string $s$, let $\overline{s}$ be its Watson-Crick reverse complement. Suppose $h:\Sigma^*\to\mathbb{Z}$ is an arbitrary string hash function. You can define $$ \tilde{h}(s)\triangleq \min\{h(s),h(\overline{s})\} $$ Then for any DNA string $s$ $$ \tilde{h}(s)=\...


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


6

Chromosomer only exists for Python 2. You should thus be able to install it via pip2 install chromosomer But: Your Python installation is a bit screwed up: there are Python 2 and Python 3, which are unfortunately incompatible. On your system python and pip seem to be aliases for Python 3, which I’d strongly recommend against (lots of tools will break). ...


6

Using SeqIO.index rather than SeqIO.parse lets you read all the records into a dict, from which you can then extract the raw entry: from Bio import SeqIO record_dict = SeqIO.index('Input.txt', 'swiss') for key in record_dict: print(record_dict.get_raw(key).decode()) Now you should be able to apply your test for a transmembrane protein to each entry, ...


6

You will have to define the header either from an existing VCF file or hardcoded into you python script. Then write the header the the output VCF file then write the dataframe to the same file with the mode options set to 'a' to append to the end of the file. You will need to order the dataframe columns to match the VCF spec and insert the empty/ default ...


6

255 is the default maximum size of a Counttable in khmer. You want to do the following: import khmer counts = khmer.Counttable(31, 1e7, 1) # Feel free to change the parameters counts.set_use_bigcount(True) The last line increases the maximum value from 255 to 65535. I don't think there's any way to currently go above that without changing the khmer code (...


6

One can get it to work by using SeqIO.InsdcIO.GenBankCdsFeatureIterator: from Bio import SeqIO file_name = 'NC_000913.3.gb' # stores all the CDS entries all_entries = [] with open(file_name, 'r') as GBFile: GBcds = SeqIO.InsdcIO.GenBankCdsFeatureIterator(GBFile) for cds in GBcds: if cds.seq is not None: cds.id = cds.name ...


6

Your setup: import pandas as pd dict1 = {0:['chr1','chr1','chr1','chr1','chr2'], 1:[1, 100, 150, 900, 1], 2:[100, 200, 500, 950, 100], 3:['feature1', 'feature2', 'feature3', 'feature4', 'feature4'], 4:[0, 0, 0, 0, 0], 5:['+','+','-','+','+']} df1 = pd.DataFrame(dict1) print(df1) ## 0 1 2 3 4 5 ## 0 chr1 1 100 ...


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

BioPython has some good tools for processing reads and alignments. http://biopython.org/DIST/docs/tutorial/Tutorial.html There is a python library wrapping samtools so many of the samtools calls can be used directly as python objects and calls https://pysam.readthedocs.io/en/latest/ I would use subprocess to call the aligner and specify the output to a bam ...


5

You're reinventing bedtools intersect (or bedops), for which there's already a convenient python module: from pybedtools import BedTool s3 = BedTool('s3.bed') s4 = BedTool('s4.bed') print(s4.intersect(s3, wa=True, wb=True, F=1)) The wb=True is equivalent to -wb with bedtools intersect on the command line. Similarly, F=1 is the same as -F 1.


5

I also wrote a package to create various plots from Oxford Nanopore sequencing data and alignments: NanoPlot. It can be installed through pip (see also the README on Github). In addition to multiple plots also a limited NanoStats output is created (see also NanoStat). Data can be presented using: A fastq file (optionally compressed) A bam file The ...


5

I whipped this little script up using the KEGG API: #!/usr/bin/env python3 import urllib.request import re import sys pathway = 'hsa00010' # glycolysis url = "http://rest.kegg.jp/get/" + pathway with urllib.request.urlopen(url) as f: lines = f.read().decode('utf-8').splitlines() want = 0 for line in lines: fields = line.split() ##...


5

using the REST API of togows: one can fetch a json structure of the genes associated to a kegg pathway. eg: http://togows.org/entry/kegg-pathway/hsa00010/genes.json [ { "3101": "HK3; hexokinase 3 [KO:K00844] [EC:2.7.1.1]", "3098": "HK1; hexokinase 1 [KO:K00844] [EC:2.7.1.1]", "3099": "HK2; hexokinase 2 [KO:K00844] [EC:2.7.1.1]", "80201": "...


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