前言

新的需求

上周我那朋友又跟我说了新的需求，然而我没空。这次他的需求是希望能统计AGCT的精确准确率，通俗点意思是从序列中找出N个连续的碱基跟目标序列比对，看多少个连续碱基就测不准了。好吧，又是新的方法。

新的问题

我考，我想，你一天一个需求，我总不能天天改代码吧，之前写的代码咋整啊，然后决定，把！我！家！的！房！子！打！理！下！，是的，就是修修剪剪缝缝补补，然后做成下面这样

#coding=utf-8
import pysam
import re

class Gene:
    def __init__(self,samfile):        ##继承外部全局变量samfile
        self.samfile =samfile
    def GC_squence_read(self,number):  ##讀取每行的數據,傳入讀取行數
        n = 0
        for line in samfile:
            n += 1
            print(n,line)
            if n == number:
                break
    def GC_total_read(self,squence_number): ##用于计算GC含量的功能模块
        number = 0
        total_GC_MAX = 0 #GC含量最大值計數
        total_GC_MIN = 1 #GC含量最小值計數
        number_lose= []  #提取失敗序列統計
        for line in samfile:
            number += 1
            if number < squence_number:
                squence = (re.findall("[AGCT][AGCT][AGCT].*[AGCT][AGCT][AGCT]",str(line))) #提取目標序列
                str_squence = " ".join(squence)
                total_squence =len(str_squence) #統計序列長度
                if total_squence < 400: #這裏是因爲曬篩除則中有一些提取出問題的序列,原因不明
                    print("序列長度：",total_squence)
                    find_G = re.findall("G",str_squence)
                    total_fin_G =len(find_G) #統計G數量
                    print("C數量",total_fin_G)
                    find_C = re.findall("C",str_squence)
                    total_fin_C = len(find_C) #統計C數量
                    print("G數量",total_fin_C)
                    total_GC = (int(total_fin_G )+ int(total_fin_C))/int(total_squence) #計算GC含量
                    print(number,total_GC)
                    if total_GC > total_GC_MAX:
                        total_GC_MAX = total_GC #判断GC含量
                        number_MAX = number
                    if total_GC_MIN > total_GC: #判断GC含量最小值
                        total_GC_MIN = total_GC
                        number_MIN = number
                else:
                    number_lose.append(number)

            if number == squence_number:
                break
        print("統計條數%s,最大GC含量%s,編號%s,最小GC含量%s,編號%s,"%(squence_number,total_GC_MAX,number_MAX,total_GC_MIN,number_MIN))
        print("提取失敗序列合計%s條,編號爲%s："%(len(number_lose),number_lose))

if __name__ == '__main__':
        path_in = "/home/charmflystar/桌面/LibPrep70_rawlib.bam"#導入BAM文件地址
        samfile = pysam.AlignmentFile(path_in, "rb")
        f1 = Gene(samfile) #實例化多個BAM文件
        f1.GC_total_read(1000000)   #參數爲需要統計的條數,以统计前100W条序列为例
        #f1.GC_squence_read(20)  #參數爲讀取的條數,用於測試

写成一个类和对象，全都是现学现用的玩意，话说编程真神奇，书本没看懂，我觉得应该是这样写，结果居然就成功了，总觉得自己有的厉害。
不惯怎么说，之前的程序感觉就不像程序，现在，总算有点像

茅房

了
咦，话说今天我应该做什么来着？