数据来源：http://hgdownload.cse.ucsc.edu/goldenPath/hg19/bigZips/chromFa.tar.gz 下载.gz数据后解压
    思路：用re模块中的正则匹配和find_all

    1. ## 目标：输出每条染色提的长度，gc含量，每一个碱基的个数
    2. import time
    3. import os
    4. import re
    6. Dir = r'/Users/Python/training/weekly/chromFa/'
    8. def job(file):
    9.     with open(Dir+file) as ch:
    10.         global A_content,G_content,C_content,T_content,N_content,result
    11.         result = {}
    12.         A_content = 0
    13.         G_content = 0
    14.         C_content = 0
    15.         T_content = 0
    16.         N_content = 0
    17.         find_A = re.compile('A',re.I) ## 忽略大小写
    18.         find_G = re.compile('G',re.I)
    19.         find_C = re.compile('C', re.I)
    20.         find_T = re.compile('T', re.I)
    21.         find_N = re.compile('N', re.I)
    22.         for line in ch:
    23.             if line.startswith(">"):
    24.                 pass
    25.             # print(line)
    26.             # print(len(line))
    27.             A = len(re.findall(find_A,line)) ## 找到每一行A碱基的个数，逐个相加
    28.             A_content += A
    29.             G = len(re.findall(find_G, line))
    30.             G_content += G
    31.             C = len(re.findall(find_C, line))
    32.             C_content += C
    33.             T = len(re.findall(find_T, line))
    34.             T_content += T
    35.             N = len(re.findall(find_N, line))
    36.             N_content += N
    37.             result['A_content'] = A_content
    38.             result['G_content'] = G_content
    39.             result['C_content'] = C_content
    40.             result['T_content'] = T_content
    41.             result['N_content'] = N_content
    42.         return (result)
    44. ## 用ch1试一下，返回一个dic，key是碱基种类，value是碱基个数，用时约70s
    46. # ch1 = job('chr1.fa')
    47. # print(ch1)
    48. # {'A_content': 65570891, 'G_content': 47016562, 'C_content': 47024413, 'T_content': 65668756, 'N_content': 23970000}
    51. ## 准备所有要计算的染色体
    52. files = []
    53. for i in (list(range(1,23))+['M','X','Y']):
    54.     file = 'chr' +str(i) + '.fa'
    55.     files.append(file)
    57. ## 4个核，约2min多
    58. from concurrent.futures import ProcessPoolExecutor
    59. with ProcessPoolExecutor(max_workers=4) as pool:
    60.     content = pool.map(job,files)
    61.     results = list(zip(files,content))
    62.     print(results)
    64. ## 计算
    65. res = []
    66. for i in range(len(results)):
    67.     chrom = results[i][0].split(".")[0]
    68.     dic = results[i][1]
    69.     dic['length'] = sum(dic.values())
    70.     dic['GC_percent'] = round((dic['G_content'] +dic['C_content'])/(dic['length']-dic['N_content']),2)
    71.     dic['N_percent'] = round(dic['N_content']/dic['length'],2)
    72.     dic['chromosome'] =chrom
    73.     res.append(dic)
    74. #print(res)
    76. ## 因为是dict，很容易放入pandas的dataframe中
    77. import pandas as pd
    78. df = pd.DataFrame.from_dict(res)

    Screen Shot 2021-11-28 at 11.16.56 PM.png
    然后可视化一下

    1. import matplotlib.pyplot as plt
    2. n=25
    3. X = np.arange(n)
    4. Y1 = df['GC_percent']*100
    5. Y2 = df['length']/1000000
    6. name = df['chromosome']
    7. name = [x[3::] for x in name]
    8. print(name)
    10. plt.bar(X,-Y1,facecolor='#9999ff',edgecolor='white')
    11. plt.bar(X,Y2,facecolor='#ff9999',edgecolor='white')
    12. Y1 = [int(x) for x in Y1]
    13. Y2 = [int(x) for x in Y2]
    15. print(Y1)
    16. print(Y2)
    17. for x,y in zip(X,Y2):
    18.     plt.text(x+0.1,y+0.05,'%d'%y,ha='center',va='bottom',fontsize=8)
    19. for x, y in zip(X, Y1):
    20.     plt.text(x + 0.1, -y - 0.08, '%d' %y, ha='center', va='top',fontsize=8)
    22. for x,y in zip(X,name):
    23.     plt.text(x+0.1,0,'%s'%y,ha='center',va='bottom',fontsize=10)
    24. plt.xlim(-.5,n)
    25. plt.ylim(-100,300)
    26. plt.title('choromsome length and GC percent')
    27. plt.yticks([-100,-50,0,100,200,300],
    28.            [r'100%',r'50%',r'0',r'100',r'200',r'300'])
    29. plt.ylabel('  percent & length(M)')
    30. plt.show()

    Screen Shot 2021-11-28 at 11.18.45 PM.png
    还有更快的方法，再试一下～