且构网

分享程序员开发的那些事...
且构网 - 分享程序员编程开发的那些事

Rosalind工具库:使用Biopython处理生物数据

更新时间:2022-02-05 12:30:34

DNA互补链 Complementing a Strand of DNA

根据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)

氨基酸翻译 Protein Translation

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

Rosalind工具库:使用Biopython处理生物数据
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

FASTQ质量分布 Base Quality Distribution

对于测序仪器而言,一般刚开始的几个碱基由于机器刚启动会存在不稳定导致质量不够,后面几个碱基会因为反应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)