python的生物学库_Rosalind工具库：使用Biopython处理生物数据

根据Waston和Crick的双螺旋学说，DNA是存在两条链，并且根据A-T, C-G反向互补。

由于DNA存在两条链，因此在分析基因组DNA序列时，我们需要要对这两条链分别处理。处理这个任务的工具有很多种，比如说EMBOSS/revseq.

问题：给定FASTA格式的DNA序列，判断在这些序列中，有多少条序列的反向互补链和原来的链是相同的。

解答：本来想用命令行工具，但是仔细一看这个题目居然还由条件判断部分，那就只能上biopython了

from Bio import SeqIO

cout = 0

for record in SeqIO.parse("rosalind_rvco.txt","fasta"):

if str(record.seq) == str(record.reverse_complement().seq):

count += 1

print(count)

给定一条DNA序列，我们希望知道这条序列是否是基因组的编码区，如果是编码区，那么我们就需要知道它会翻译成什么样的氨基酸。

DNA-AA

尽管大部分生物享有相同的翻译机制，比如说起始密码子为(AUG)，终止密码子为(UAA,UAG,UGA), 但在部分物种和含DNA的细胞器上存在着遗传密码变体(genetic code variation)，比如说脊椎动物的线粒体将AGA和AGG当作终止密码子。

问题：给定一条DNA序列和AA序列，判断从DNA到AA的翻译使用了哪个翻译表？

这个题目非常简单，但是却让我发现原来从DNA到AA的翻译有那么多方式，NCBI上的编号显示一共有将近30种。解题策略就是用Biopython使用不同的密码表去翻译，然后判断结果和AA是否相同。

from Bio.Seq import translate

from Bio.Seq import CodonTable

def check_table(dna,aa):

for i in CodonTable.ambiguous_dna_by_id:

if translate(dna,table=i,to_stop=True) == aa:

return i

data = [line.strip() for line in open("rosalind_ptra.txt")]

check_table(data[0],data[1])

序列质量分布

这题是关于高通量测序下机数据（FASTQ）的质量分布，题目非常的简单，就是对给定的序列统计低于某个阈值的read数，原本打算使用fastp，但是做了好几次都出错。后来发现发现是fastq的phred质量编码出了问题，在历史上曾经出现过三种编码方式

fastq-sanger: Phread+33

fastq-solexa: Phred+64

fastq-illumina: Phred+64

考虑到Rosalind这道题目的还是熟悉Biopython，为了通过就只能用Biopython了

from Bio import SeqIO

for record in SeqIO.parse("rosalind_phre.txt","fastq-sanger"):

total = sum(record.letter_annotations['phred_quality'])

mean = total / len(record.letter_annotations['phred_quality'])

if mean < 27:

count += 1

对于测序仪器而言，一般刚开始的几个碱基由于机器刚启动会存在不稳定导致质量不够，后面几个碱基会因为反应DNA聚合酶活性不够等因素导致质量不佳。

FastQC质控过程就有一个指标判断每个位置上的碱基分布是否符合要求，称之为"Per Base Sequence Quality"

问题： 统计fastq文件每个位置上的碱基的质量平均值，要求给出低于平均值的位置。

答案： 使用biopython读取序列，使用numpy进行矩阵运算

from Bio import SeqIO

import numpy as np

lists = []

for record in SeqIO.parse("rosalind_bphr.fq","fastq-sanger")

lists.append(record.letter_annotations['phred_quality'])

means = np.mean(np.array(lists), axis=0)

np.sum(means < 24)