【Python】生信编程实战

作者: caokai001 | 来源:发表于2019-01-21 02:00 被阅读29次

第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

相关文章

网友评论

    本文标题:【Python】生信编程实战

    本文链接:https://www.haomeiwen.com/subject/frlqjqtx.html