前两节我们已经学习了利用Biostrings包进行基因序列的基本操作,模板匹配,搜索回文结构。接着上一节的内容我们今天来学习应用Biostrings做生物信息数据的序列对比。
上一节的内容:应用Biostrings处理生物信息数据——基础篇(二)
本期的代码:
1、序列对比
# 用AAString函数生成一个"AAString"对象aa1
aa1# 用AAString函数生成一个"AAString"对象aa2
aa2#序列全局比对,应用矩阵"BLOSUM62"打分,要求gapopen罚分为3。
needwunsQS(aa1, aa2, "BLOSUM62", gappen =3)#用DNAString函数生成一个"DNAString"对象dnal。
dnal#用DNAString函数生成-一个"DNAString"对象dna2。
dna2#构建4X4的矩阵DNA打分矩阵。
mat#序列全局比对,应用矩阵mat打分,要求gapopen罚分为0
needwunsQS(dnal, dna2, mat, gappen =0)2、读写序列文件(Fasta 和Fastq格式)
#指定文件的目录(Biostrings安装目录中的extdata子目录)和文件名(someORFfa)。
filepath#读取FASTA文件。
x#查看FASTA文件的内容。
x#命名输出文件
out1#把序列输出到文件out1,格式还是FASTA。
writeXStringSet(x, out1)#指定文件的目录和文件名(s_1_sequence.txt), 这个是FASTQ文件。
filepath#显示上面FASTQ文件中的数据信息。
fastq.geometry(filepath)#表示数据有256条reads, 测序长度是36bp#读取FASTQ文件。
x#查看FASTQ文件,结果不在这里显示。
x#从第1号染色体上按照固定长度(50bp)依次取短序列(read) 的起点(向量)。
sw_start#从起点开始取read,每个长度为10bp, 注意sw的格式是"XStringViews"。
sw#变量sw的格式从"XStringViews"转换为"XStringSet"。
my_fake_shorreads#按照"ID"加开头6个数字的格式,得到一组新名称。
my_fake_ids#用新名称替换旧名称。
names(my_fake_shorreads)#查看第50000到50005条数据,结果不在这里显示。
my_fake_shorreads[500000:500005]#命名输出文件。
out2#把序列输出到文件out2,格式是FASTQ,但是缺少质量信息。
writeXStringSet(my_fake_shorreads, out2, format ="fastq")#产生质量信息。
my_fake_quals =#查看my_fake_quals 内容。
my_fake_quals#命名输出文件。
out3#把序列输出到文件out3,格式还是FASTQ,这次 含有质量信息。
writeXStringSet(my_fake_shorreads, out3, format ="fastq", qualities = my_fake_quals)