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序列,我们希望知道这条序列是否是基因组的编码区,如果是编码区,那么我们就需要知道它会翻译成什么样的氨基酸。
尽管大部分生物享有相同的翻译机制,比如说起始密码子为(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)