第0题- 对fasta文件取互补序列
知识要点
1:对字典对象的学习;
2:replace方法一次只能完成一种替换,即某字符串替换成某字符串;(并不像perl中的tr/ATCG/TAGC/,能够一次完成多种替换)
① 使用字典替换
##cat test2.fa
>NM_001011874 gene=Xkr4 CDS=151-2091
gcggcggcgggcgagcgggcgctggagtaggagctggggagcggcgcggccggggaaggaagccagggcg
aggcgaggaggtggcgggaggaggagacagcagggacaggTGTCAGATAAAGGAGTGCTCTCCTCCGCTG
CCGAGGCATCATGGCCGCTAAGTCAGACGGGAGGCTGAAGATGAAGAAGAGCAGCGACGTGGCGTTCACC
CCGCTGCAGAACTCGGACAATTCGGGCTCTGTGCAAGGACTGGCTCCAGGCTTGCCGTCGGGGTCCGGAG
>NM_001195662 gene=Rp1 CDS=55-909
AAGCTCAGCCTTTGCTCAGATTCTCCTCTTGATGAAACAAAGGGATTTCTGCACATGCTTGAGAAATTGC
AGGTCTCACCCAAAATGAGTGACACACCTTCTACTAGTTTCTCCATGATTCATCTGACTTCTGAAGGTCA
AGTTCCTTCCCCTCGCCATTCAAATATCACTCATCCTGTAGTGGCTAAACGCATCAGTTTCTATAAGAGT
GGAGACCCACAGTTTGGCGGCGTTCGGGTGGTGGTCAACCCTCGTTCCTTTAAGACTTTTGACGCTCTGC
TGGACAGTTTATCCAGGAAGGTACCCCTGCCCTTTGGGGTAAGGAACATCAGCACGCCCCGTGGACGACA
CAGCATCACCAGGCTGGAGGAGCTAGAGGACGGCAAGTCTTATGTGTGCTCCCACAATAAGAAGGTGCTG
# python3代码
complement = {'A': 'T', 'G': 'C', 'C': 'G', 'T': 'A', 'N': 'N'}
with open("test2.fa") as f:
for line in f:
line_strip=line.strip()
if line_strip.startswith(">"):
print("\n"+line_strip,end="\n")
else:
seq=""
for i in line_strip.upper():
seq=seq+complement[i]
print(seq,end="")
#运行结果
>NM_001011874 gene=Xkr4 CDS=151-2091
CGCCGCCGCCCGCTCGCCCGCGACCTCATCCTCGACCCCTCGCCGCGCCGGCCCCTTCCTTCGGTCCCGCTCCGCTCCTCCACCGCCCTCCTCCTCTGTCGTCCCTGTCCACAGTCTATTTCCTCACGAGAGGAGGCGACGGCTCCGTAGTACCGGCGATTCAGTCTGCCCTCCGACTTCTACTTCTTCTCGTCGCTGCACCGCAAGTGGGGCGACGTCTTGAGCCTGTTAAGCCCGAGACACGTTCCTGACCGAGGTCCGAACGGCAGCCCCAGGCCTC
>NM_001195662 gene=Rp1 CDS=55-909
TTCGAGTCGGAAACGAGTCTAAGAGGAGAACTACTTTGTTTCCCTAAAGACGTGTACGAACTCTTTAACGTCCAGAGTGGGTTTTACTCACTGTGTGGAAGATGATCAAAGAGGTACTAAGTAGACTGAAGACTTCCAGTTCAAGGAAGGGGAGCGGTAAGTTTATAGTGAGTAGGACATCACCGATTTGCGTAGTCAAAGATATTCTCACCTCTGGGTGTCAAACCGCCGCAAGCCCACCACCAGTTGGGAGCAAGGAAATTCTGAAAACTGCGAGACGACCTGTCAAATAGGTCCTTCCATGGGGACGGGAAACCCCATTCCTTGTAGTCGTGCGGGGCACCTGCTGTGTCGTAGTGGTCCGACCTCCTCGATCTCCTGCCGTTCAGAATACACACGAGGGTGTTATTCTTCCACGAC
② 使用sed 来实现快速替换
sed 's/[aA]/T/gI' file
③ 使用Biopython来获取反向序列
对于核苷酸序列,你可以使用 Seq 对象内置的方法很容易地获得 Seq 的互补或反向互补序列。
from Bio.Seq import Seq
from Bio.Alphabet import IUPAC
my_seq = Seq("GATCGATGGGCCTATATAGGATCGAAAATCGC", IUPAC.unambiguous_dna)
#my_seq
Seq('GATCGATGGGCCTATATAGGATCGAAAATCGC', IUPACUnambiguousDNA())
#my_seq.complement()
Seq('CTAGCTACCCGGATATATCCTAGCTTTTAGCG', IUPACUnambiguousDNA())
#my_seq.reverse_complement()
Seq('GCGATTTTCGATCCTATATAGGCCCATCGATC', IUPACUnambiguousDNA())
在前面的方法中,使用切片的-1的步长可以很容易的获取一个 Seq 对象的反向序列。
my_seq[::-1]
Seq('CGCTAAAAGCTAGGATATATCCGGGTAGCTAG', IUPACUnambiguousDNA())
第1题- 序列长度分布统计
对任意fastq文件进行序列长度分布统计:总长、平均值、最长、最短长度分别是多少?
**知识点**:
列表、文件读入、打印输出、导入包并使用基函数、
## cat test1.fq
@SRR400264.2496 HWI-ST216_0180:3:1101:11339:2337 length=36
CTGCCCCCGCTAACCGGCTTTTTGCCCAAATGGGCC
+SRR400264.2496 HWI-ST216_0180:3:1101:11339:2337 length=36
HHHHHHHHHHHHHHHHHHGHHHHGHHGHHGHHFEHH
@SRR400264.2497 HWI-ST216_0180:3:1101:11298:2341 length=35
AGCTTTTTTTTTTCTTTTTCTTTTTTGAGATGGCA
+SRR400264.2497 HWI-ST216_0180:3:1101:11298:2341 length=35
HHHHHHHHHHHHHGHHHHHHHHHHHHDFHAF=GGF
@SRR400264.2498 HWI-ST216_0180:3:1101:11256:2344 length=33
CTGTCTTCTTCCCCAGTTCATCTGCATCTCGTT
+SRR400264.2498 HWI-ST216_0180:3:1101:11256:2344 length=33
HHHHHHHHHHHHHHHHHGGHGHHGHGHHFHHHH
@SRR400264.2499 HWI-ST216_0180:3:1101:11384:2357 length=30
ATCGATAACTAAACTTCCTTTCCATATGAC
+SRR400264.2499 HWI-ST216_0180:3:1101:11384:2357 length=30
HHHHHHHHHHHHHHHHHHHHHHHHHHHHHH
@SRR400264.2500 HWI-ST216_0180:3:1101:11352:2366 length=35
CCGTAACTCTTATAAGCTAGCTTATATAAGAGCTT
+SRR400264.2500 HWI-ST216_0180:3:1101:11352:2366 length=35
HHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHGHG
代码如下
lens=[]
len_sum=0
###读入文件
with open("test2.fa") as f:
count=0
for line in f:
count+=1
line_strip=line.strip()
if(count %4)==0: ## 判断是否是质量行,接着用len()统计长度
length=len(line.strip())
lens.append(length) ## 每遇到一次质量行统计长度后把长度值加到lens列表中
#使用sum()/mean()/min()函数对lens列表进行统计:
Sum=sum(lens)
Mean= mean(lens)
Min=min(lens)
# 使用len()函数对lens列表计算元素个数
Num=len(lens)
print("Sum",Sum,sep=":")
print("Mean", Mean, sep=":")
print("Min", Min, sep=":")
print("Num", Num, sep=":")
## result:
Sum:169
Mean:33.8
Min:30
Num:5
第2题- fa/fq文本处理
##2.1写程序 splitName.py, 读入test2.fa, 并取原始序列名字第一个空格前的名字为处理后的序列名字,输出到屏幕
with open("test2.fa") as f:
for line in f:
if line.startswith(">"):
print(line.split(" ")[0])
else:
print(line.strip())
#结果
>NM_001011874
gcggcggcgggcgagcgggcgctggagtaggagctggggagcggcgcggccggggaaggaagccagggcg
aggcgaggaggtggcgggaggaggagacagcagggacaggTGTCAGATAAAGGAGTGCTCTCCTCCGCTG
CCGAGGCATCATGGCCGCTAAGTCAGACGGGAGGCTGAAGATGAAGAAGAGCAGCGACGTGGCGTTCACC
CCGCTGCAGAACTCGGACAATTCGGGCTCTGTGCAAGGACTGGCTCCAGGCTTGCCGTCGGGGTCCGGAG
>NM_001195662
AAGCTCAGCCTTTGCTCAGATTCTCCTCTTGATGAAACAAAGGGATTTCTGCACATGCTTGAGAAATTGC
AGGTCTCACCCAAAATGAGTGACACACCTTCTACTAGTTTCTCCATGATTCATCTGACTTCTGAAGGTCA
AGTTCCTTCCCCTCGCCATTCAAATATCACTCATCCTGTAGTGGCTAAACGCATCAGTTTCTATAAGAGT
GGAGACCCACAGTTTGGCGGCGTTCGGGTGGTGGTCAACCCTCGTTCCTTTAAGACTTTTGACGCTCTGC
TGGACAGTTTATCCAGGAAGGTACCCCTGCCCTTTGGGGTAAGGAACATCAGCACGCCCCGTGGACGACA
CAGCATCACCAGGCTGGAGGAGCTAGAGGACGGCAAGTCTTATGTGTGCTCCCACAATAAGAAGGTGCTG
#2.2
网友评论