我有一个文本文件,其中包含蛋白质序列(200个序列),如下所示.
>ptn1
AAGHM
>ptn2
MGLKKRR
我需要为seqence的每个字符提供以下值,并且必须找到每个序列的平均值.
A= 0.2, G= 0.5, L=0.14, M= 0.70, R= 0.55, C=0.48, H= 1.00 , K=0.4
期望的输出
ptn1 - 0.52
ptn2 - 0.462
我怎么能用awk或python做到这一点？
您的建议将不胜感激
我想从pdb文件中提取链.我有一个名为pdb.txt的文件,其中包含pdb ID,如下所示.前四个字符表示PDB ID,最后一个字符是链ID.
1B68A
1BZ4B
4FUTA
我想1)逐行读取文件2)从相应的PDB文件下载每个链的原子坐标.
3)将输出保存到文件夹.
我使用以下脚本来提取链.但是此代码仅打印来自pdb文件的A链.
for i in 1B68 1BZ4 4FUT
do
wget -c "http://www.pdb.org/pdb/download/downloadFile.do?fileFormat=pdb&compression=NO&structureId="$i -O $i.pdb
grep ATOM $i.pdb | grep 'A' > $i\_A.pdb
done
我有一个文本文件,如下所示
ATOM 920 CA GLN A 203 39.292 -13.354 17.416 1.00 55.76 C
ATOM 929 CA HIS A 204 38.546 -15.963 14.792 1.00 29.53 C
ATOM 939 CA ASN A 205 39.443 -17.018 11.206 1.00 54.49 C
ATOM 947 CA GLU A 206 41.454 -13.901 10.155 1.00 26.32 C
ATOM 956 CA VAL A 207 43.664 -14.041 13.279 1.00 40.65 C
.
.
.
ATOM 963 CA GLU A 208 45.403 -17.443 13.188 1.00 40.25 C
我想计算两个α碳原子之间的距离,即计算第一和第二原子之间的距离,然后计算第二和第三原子之间的距离等等......两个原子之间的距离可以表示为:distance = sqrt((x1-x2)^2+(y1-y2)^2+(z1-z2)^2) …
我正在研究.smiles文件..smiles文件的文件结构是:http://en.wikipedia.org/wiki/Chemical_file_format#SMILES
我想从微笑文件中获取所有原子.这意味着如果存在单个"C"原子,则意味着将有4个"H"原子连接到它们.
我在搜索时发现python中有一些模块可以解析微笑格式,但它们不会给出支持的氢原子.(例如:它们只给'C'而不是其他4'H'原子连接到'C'原子)
如何使用python找到所有原子,包括连接的'H'原子.
需要转换为包含连接的'H'原子的所有原子的smiles文件示例:
[H]OC([H])([H])[C@@]1([H])C([H])=C([H])[C@@]([H])(n2c([H])nc3c(nc(nc23)N([H])[H])N([H])C2([H])C([H])([H])C2([H])[H])C1([H])[H]
先感谢您.
我开发了以下代码来计算对齐中相同站点的数量.不幸的是代码很慢,而且我必须在数百个文件上进行迭代,处理超过1000个对齐需要将近12个小时,这意味着代码的速度要快十倍.任何帮助,将不胜感激:
import os
from Bio import SeqIO
from Bio.Seq import Seq
from Bio import AlignIO
from Bio.SeqRecord import SeqRecord
from Bio.Alphabet import generic_dna
from Bio.Align import MultipleSeqAlignment
import time
a = SeqRecord(Seq("CCAAGCTGAATCAGCTGGCGGAGTCACTGAAACTGGAGCACCAGTTCCTAAGAGTTCCTTTCGAGCACTACAAGAAGACGATTCGCGCGAACCACCGCAT", generic_dna), id="Alpha")
b = SeqRecord(Seq("CGAAGCTGACTCAGTGGGCGGAGTCACTGAAACTGGAGCACCAGTTCCTCAGAGTCCCCTTCGAGCACTACAAGAAGACAATTCGTGCGAACCACCGCAT", generic_dna), id="Beta")
c = SeqRecord(Seq("CGAAGCTGACTCAGTTGGCAGAATCACTGAAACTGGAGCACCAGTTCCTCAGAGTCCCCTTCGAGCACTACAAGAAGACGATTCGTGCGAACCACCGCAT", generic_dna), id="Gamma")
d = SeqRecord(Seq("CGAAGCTGACTCAGTTGGCAGAGTCACTGAAACTGGAGCACCAGTTCCTCAGAGTCCCCTTCGAGCACTACAAGAAGACGATTCGTGCGAACCACCGCAT", generic_dna), id="Delta")
e = SeqRecord(Seq("CGAAGCTGACTCAGTTGGCGGAGTCACTGAAACTGGAGCACCAGTTCCTCAGAGTCCCCTTCGAGCACTACAAGAAGACGATTCGTGCGAACCACCGCAT", generic_dna), id="Epsilon")
align = MultipleSeqAlignment([a, b, c], annotations={"tool": "demo"})
start_time = time.time()
if len(align) != 1:
for n in range(0,len(align[0])):
n=0
i=0
while n<len(align[0]): #part that needs to …我有一个fasta文件(下面提到的第一个序列)有很长的描述.我需要选择具体的描述字段.当我使用下面的代码; 整个描述进入字符串.
from Bio import SeqIO
for record in SeqIO.parse("geneTemp.fasta", "fasta") :
id=record.id
desc=record.description
print desc
有没有简单的方法将描述字段(使用biopython库)放入数组并选择特定字段而不将描述转换为字符串并吐出字符串？
代码输出
Python 2.7 (r27:82500, Sep 16 2010, 18:03:06)
[GCC 4.5.1 20100907 (Red Hat 4.5.1-3)] on localhost.localdomain, Standard
>>> FBgn0197520 type=gene; loc=scaffold_12855:complement(6241650..6242111); ID=FBgn0197520; name=Dvir\GJ10233; dbxref=FlyBase_Annotation_IDs:GJ10233,FlyBase:FBgn0197520,GLEANR:dvir_GLEANR_10171,EntrezGene:6632532,GB_protein:EDW59542,FlyMine:FBgn0197520,OrthoDB4.Arthropods:FBgn0242841,OrthoDB4.Arthropods:FBgn0213090,OrthoDB4.Arthropods:FBgn0190974,OrthoDB4.Arthropods:FBgn0165423,OrthoDB4.Arthropods:FBgn0247590,OrthoDB4.Arthropods:FBgn0149779,OrthoDB4.Arthropods:FBgn0146205,OrthoDB4.Arthropods:FBgn0017456,OrthoDB4.Arthropods:FBgn0126736,OrthoDB4.Arthropods:FBgn0117264,OrthoDB4.Arthropods:FBgn0094317; MD5=0b7e859d2a6eca028ffd16b964835705; length=462; release=r1.2; species=Dvir;
loc=scaffold_12855:complement(6241650..6242111)
来自fasta文件的序列之一.
>FBgn0207418 type=gene; loc=scaffold_12875:complement(14361770..14363857); ID=FBgn0207418; name=Dvir\GJ20278; dbxref=FlyBase_Annotation_IDs:GJ20278,FlyBase:FBgn0207418,GLEANR:dvir_GLEANR_5721,EntrezGene:6625684,GB_protein:EDW61510,FlyMine:FBgn0207418,OrthoDB4.Arthropods:NV16422,OrthoDB4.Arthropods:LH16819,OrthoDB4.Arthropods:ISCW000548,OrthoDB4.Arthropods:FBgn0239668,OrthoDB4.Arthropods:FBgn0219970,OrthoDB4.Arthropods:FBgn0181866,OrthoDB4.Arthropods:FBgn0175499,OrthoDB4.Arthropods:FBgn0080765,OrthoDB4.Arthropods:FBgn0155230,OrthoDB4.Arthropods:FBgn0141947,OrthoDB4.Arthropods:FBgn0033392,OrthoDB4.Arthropods:FBgn0127494,OrthoDB4.Arthropods:FBgn0102879,OrthoDB4.Arthropods:FBgn0090125,OrthoDB4.Arthropods:CPIJ005729,OrthoDB4.Arthropods:GB15324,OrthoDB4.Arthropods:AGAP012336,OrthoDB4.Arthropods:AAEL007395,OrthoDB4.Arthropods:PB24927,OrthoDB4.Arthropods:PHUM365660,OrthoDB4.Arthropods:GLEAN_06039; MD5=4c62b751ec045ac93306ce7c08d254f9; length=2088; release=r1.2; species=Dvir;
ATGCGTCTGCGACGCCGCTGGCATCGGCGGATGCGGCGTACAATTGAGAA
AATCTATCGCCTTAAAATGCAATCGCGCCGCAAGTTGGTTTACTTAGCCG
TATTTGGAGCACTATGCGTAATATTCTGGCTGGCTGGACAGCAGTTGCTG
ACGACTTCGAATGGTCACTACAGTAGCTACTACGGCGAAACGCATTGTGC
GCCCATTGATGCCGTATACACCTGGGTAAATGGTTCGGATCCGGATTTTA
TTGAGTCCATTAGACGCTACGATGCCAGCTACGATCCGTCGCGCTTCGAC
有人可以向我展示如何计算 DNA 序列中最长开放阅读框 (ORF) 的简单解决方案吗？ ATG是起始密码子（即，一个ORF的开始）和TAG，TGA以及TAA是终止密码子（即，一个ORF的结束）。
这是一些产生错误的代码（并使用名为 BioPython 的外部模块）：
import sys
from Bio import SeqIO
currentCid = ''
buffer = []
for record in SeqIO.parse(open(sys.argv[1]),"fasta"):
cid = str(record.description).split('.')[0][1:]
if currentCid == '':
currentCid = cid
else:
if cid != currentCid:
buffer.sort(key = lambda x : len(x[1]))
print '>' + buffer[-1][0]
print buffer[-1][1]
currentCid = cid
buffer = [(str(record.description),str(record.seq))]
else:
buffer.append((str(record.description),str(record.seq)))
buffer.sort(key = lambda x : len(x[1]))
print '>' + buffer[-1][0]
print buffer[-1][1]
是否有可能用最少的外部依赖来编写这个过程（或者至少让上面的代码工作）？ …
我有一个纽维克格式的系统发育树。我想根据终端节点的标签（因此基于物种列表）拉出一棵子树。我正在使用的树的副本可以在这里找到： http: //hgdownload.soe.ucsc.edu/goldenPath/dm6/multiz27way/dm6.27way.nh
目前我已经使用 BioPython 阅读了树，如下所示：
from Bio import Phylo
#read in Phylogenetic Tree
tree = Phylo.read('dm6.27way.nh', 'newick')
#list of species of interest
species_list = ['dm6', 'droSim1', 'droSec1', 'droYak3', 'droEre2', 'droBia2', 'droSuz1', 'droAna3', 'droBip2', 'droEug2', 'droEle2', 'droKik2', 'droTak2', 'droRho2', 'droFic2']
我如何仅提取species_list中物种的子树？
我是 Biopython 的新手（以及一般编码），我正在尝试编写一种方法，将一系列 DNA 序列（超过 80 个）翻译成蛋白质序列，在一个单独的 FASTA 文件中。我还想在正确的阅读框中找到序列。
这是我到目前为止所拥有的：
from Bio import SeqIO
from Bio.SeqRecord import SeqRecord
for record in SeqIO.parse("dnaseq.fasta", "fasta"):
protein_id = record.id
protein1 = record.seq.translate(to_stop=True)
protein2 = record.seq[1:].translate(to_stop=True)
protein3 = record.seq[2:].translate(to_stop=True)
if len(protein1) > len(protein2) and len(protein1) > len(protein3):
protein = protein1
elif len(protein2) > len(protein1) and len(protein2) > len(protein3):
protein = protein2
else:
protein = protein3
def prot_record(record):
return SeqRecord(seq = protein, \
id = ">" + protein_id, \
description = "translated sequence")
records = …我正在尝试使用已发布的所有序列来构建数据库细菌类型，以使用 bowtie2 进行映射来计算我对这个数据库的读取覆盖率，为此，我将我从 ncbi 下载的所有基因组序列合并到一个 fasta_library 中（我合并了 74 个文件在 fasta 文件中），问题是在这个 fasta 文件（我创建的库）中我有很多重复的序列，这在很大程度上影响了覆盖率，所以我问是否有任何方法可以消除重复我的 Library_File 中有，或者是否有任何方法可以在没有重复的情况下合并序列，或者是否有任何其他方法可以计算我的读取对参考序列的覆盖率
我希望我足够清楚，如果有什么不清楚的请告诉我。