今天要连发两篇文章，那是因为...算了不矫情了
下面献上我最后的作品吧

#coding=utf-8
import pysam
import re
import os

class Gene:
    def __init__(self,samfile,path_in):        ##继承外部全局变量samfile
        self.samfile = samfile
        self.path_in = path_in
    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))
    def pysam_save(self,squence_number):  #用于提取序列的功能模块,输出TXT
        number =0
        compile = re.compile("/桌面/(.*).bam")
        filename = re.findall(compile,path_in)
        print(filename)
        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)
                if len(str_squence) < 400:
                    with open("bam_squence", "a") as pysam_txt:
                        pysam_txt.write(str_squence)
                        pysam_txt.write("\n")
                        #print(os.path.abspath("bam_squence.txt"))
            else:
                ospath = "/home/charmflystar/Desktop/save_file_for_pycharm/gene_squence/bam_squence"
                ospath_new = "/home/charmflystar/Desktop/save_file_for_pycharm/gene_squence/%s"% filename
                os.rename(ospath, ospath_new)
                print("执行完毕")
                break
    def AGCT_accuracy(self,squence_number): #用于统计A\G\C\T的碱基精准度,智能化输出到TXT
        number = 0
        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)#需要转化为str
                total_squence =len(str_squence) #統計序列長度
                if total_squence < 400: #這裏是因爲曬篩除則中有一些提取出問題的序列,原因不明
                    compile1 = re.compile("[A][A][A][A][A*]") #根据需求抓取含A的连续序列
                    A_squence = str(re.findall(compile1,str_squence))#同样需转化为str,以便以文本格式保存
                    compile2 = re.compile("/桌面/(.*).bam")#1.这段代码是用来智能化命名文本名称的
                    name = re.findall(compile2, path_in)#2.这段代码是用来智能化命名文本名称的
                    filename =" ".join(name)#3.这段代码是用来智能化命名文本名称的
                    #print(A_squence)
                    if "A" in A_squence:
                        with open("bam_squence", "a") as pysam_txt: #把筛出的符合序列写到记事本中
                            pysam_txt.write(str_squence)
                            pysam_txt.write("\n")
            else:
                ospath = "/home/charmflystar/Desktop/save_file_for_pycharm/gene_squence/bam_squence"
                ospath_new = "/home/charmflystar/Desktop/save_file_for_pycharm/gene_squence/%s" % filename #4.这段代码是用来智能化命名文本名称的
                os.rename(ospath, ospath_new)
                print("执行完毕")
                break
if __name__ == '__main__':
        path_in = "/home/charmflystar/桌面/LibPrep71_rawlib.bam"  # 導入BAM文件地址
        samfile = pysam.AlignmentFile(path_in, "rb")
        f1 = Gene(samfile,path_in) #實例化多個BAM文件
        #f1.GC_total_read(1000000)   #參數爲需要統計的條數,以统计前100W条序列为例
        #f1.GC_squence_read(20)  #參數爲讀取的條數,用於測試
        #f1.pysam_save(10)    #参数为需要保存的行数,用于抽取BAM文件中的所有序列,例如,抽取前10条序列,输出TXT文本.没事别用这个,试了一下我输出200W条序列,容量达到了惊人的300MB
        #f1.AGCT_accuracy(1000000)   #参数为需要统计的条数，以统计100W条序列为例

面向百度和书本编程