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 ...
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)
...
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
No, there is currently no GFF support in biopython.
However, you can read in GFF files into python using this package, gffutils. There are also a few other packages to read/write GFF files, like gff3.
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
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 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 ...
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
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
...
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
With a k-mer size of 28 it shouldn't be finding that many matches. And the prokka results are suspicious as well. Maybe you have multiple contigs (none larger than 100kb) in that file? What is the result of
grep ^'>' fasta_file | wc -l
? This would show how many contigs you have in the file.
6
This can be done by using the "Search details" as a search term in Entrez.esearch:
import os
OPJ = os.path.join
base_dir = os.getcwd()
from Bio import Entrez
Entrez.email = "my_email_here"
query = "ANOS1[gene name] AND refseq[filter]"
handle = Entrez.esearch(term=query, db="nuccore", retmax=400)
ids = Entrez.read(handle)["IdList"]
with open(OPJ(base_dir, "...
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
You can specify send to stdout using out='-' in the Biopython wrapper.
from Bio.Blast.Applications import NcbiblastnCommandline
import pandas as pd
cline = NcbiblastnCommandline(query='seq.fna', subject='seq2.fna', outfmt=6, out='-')
output = cline()[0].strip()
rows = [line.split() for line in output.splitlines()]
cols = ['qseqid', 'sseqid', 'pident', '...
6
No PDB
Getting a crystal structure is hard work for crystallographers, even with a high throughput systems, so they make constructs with only known domains or regions of particular interest. What happened here is that only two parts have been cloned, expressed and crystallised. The rest may be crystallisable or solvable by cryoEM, but may be too disordered.
...
5
I've written a handful of programs from scratch to simulate mutations and variations in real or simulated sequences.
The trick has always been to sort the variants by genomic coordinate, apply the variant with the largest coordinate first, then apply the variant with the second largest coordinate, all the way down to the variant with the smallest coordinate....
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
You'll first need to determine the appropriate line number, which you can do with grep -m1 -nw "@31027" foo.fastq. After that, note that you can provide a line number to sed:
sed -i '123456s/CAAAAGTCACTCA/CAAAAGTCACTTA/' foo.fastq
That will do the replacement only on line 123456 and edit the file in place (the -i option).
5
You're using version 1.68 or older. Mauve support was added in 1.70.
5
Here's one way with pure Python:
#!/usr/bin/env python3
import sys
import re
with open (sys.argv[1] , 'r') as fo:
for line in fo:
fields = line.split()
if len(fields)==2:
print(fields[1],len(fields[1].replace('-','')))
Running this on your file gives:
$ foo.py file.aln
------------------------------------------------------...
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.
4
Found what I was looking for. In:
searchResultHandle = Entrez.esearch(db="protein", term="terminase large", retmax=1000)
I've added:
searchterm = "(terminase large subunit AND viruses[Organism]) AND Caudovirales AND refseq[Filter]"
searchResultHandle = Entrez.esearch(db="protein", term=searchterm, retmax=6000)
which norrowed down my searches to the ...
4
You can use finditer from python's re module. The advantage of this approach is it allows for getting the indices of those matches, which could be handy down the track.
>>> import re
>>> matches = re.finditer(r'(?=(AA))', 'AAAA')
>>> indices = [match.span(1) for match in matches]
>>> indices
[(0, 2), (1, 3), (2, 4)]
>&...
4
It's important to always consider read length and quality jointly with high-error read data, and current long-read technologies (e.g., MinION and PacBio) have high error rates. Considering read length and quality jointly will help you determine how successful the run was, how many reads were 'high quality', whether the longer reads are 'real' (or just pore ...
4
If you look for a program which would randomly introduce SNPs + short indels and then would save everything into a VCF file, DWGsim or Mason Variator could be a good choice. Then you can create a corresponding Chain file using bcftools consensus -c and transform various formats between these two coordinate systems using CrossMap.
4
Currently you could use either but a major question is which platform will others be using in the future. AFAIK Perl is only superior to Python for regex.
Based on the trend I see for new programmers and new software being released: Perl is on the way out and Python is still growing.
https://trends.google.com/trends/explore?date=all&q=bioperl,biopython
4
This usually comes down to religious issues, so let me try and steer it back to more objective grounds:
What language do you know (better)? Use the library for that one.
If you know neither and will be learning a language to use the library, the majority opinion would be that Python is easier to learn.
However, some people say that they "click" with Perl ...
4
Regading the perl vs python discussion, there is no final answer which language is better, but I have some advice for you:
Learn the language your colleagues or your advisor use. This way you are able to discuss your code with them and also get help if you run into problems.
4
Here is a sample of code that should get you started.
I created a separate get_abstracts function. Then that function get called repetitively until less than 100K ids are returned. If you wanted to get fancy this whole process could be called as a recursive function but that seemed overly complex.
BUFFER_SIZE = 100000
def get_abstracts(start, term):
...
Only top voted, non community-wiki answers of a minimum length are eligible