说点闲话
起因是因为有个朋友问我能不能处理BAM文件，我找了一下python的确有处理BAM文件的模块，在没有搞清楚能不能做的情况下，我说可以没问题，这B可装上了天。
然后，这B是哭着也要装下去了。
首先了解一下pysam和BAM是干啥的：
python读写处理模块pysam是在开发基因组相关流程或工具时，经常需要读取、处理和创建bam、vcf、bcf文件。它打包了htslib-1.3、samtools-1.3 和 bcftools-1.3的核心功能，能在编程时非常灵活的处理bam和bcf文件。
简单来说就是一个存取数据库一样的东西，利用pysam功能可以实现提取文本，OK。
一、安装环境部署
打开cmd然后
pip install pysam
C:\Users\hasee>pip install pysam
Collecting pysam
Using cached https://files.pythonhosted.org/packages/73/59/c319f1bde3019bbce4583cecb12b9e3e52ffcfbe2c96d8b1fb131c0d4fb7/pysam-0.15.1.tar.gz
Complete output from command python setup.py egg_info:
'.' 不是内部或外部命令，也不是可运行的程序
或批处理文件。
'.' 不是内部或外部命令，也不是可运行的程序
或批处理文件。
# pysam: no cython available - using pre-compiled C
# pysam: htslib mode is shared
# pysam: HTSLIB_CONFIGURE_OPTIONS=None
# pysam: htslib configure options: None
Traceback (most recent call last):
File "<string>", line 1, in <module>
File "C:\Users\hasee\AppData\Local\Temp\pip-install-pgf7vjf2\pysam\setup.py", line 223, in <module>
htslib_make_options = run_make_print_config()
File "C:\Users\hasee\AppData\Local\Temp\pip-install-pgf7vjf2\pysam\setup.py", line 69, in run_make_print_config
stdout = subprocess.check_output(["make", "-s", "print-config"])
File "c:\python37\lib\subprocess.py", line 376, in check_output
**kwargs).stdout
File "c:\python37\lib\subprocess.py", line 453, in run
with Popen(*popenargs, **kwargs) as process:
File "c:\python37\lib\subprocess.py", line 756, in __init__
restore_signals, start_new_session)
File "c:\python37\lib\subprocess.py", line 1155, in _execute_child
startupinfo)
FileNotFoundError: [WinError 2] 系统找不到指定的文件。
----------------------------------------
Command "python setup.py egg_info" failed with error code 1 in C:\Users\hasee\AppData\Local\Temp\pip-install-pgf7vjf2\pysam\
You are using pip version 18.0, however version 18.1 is available.
You should consider upgrading via the 'python -m pip install --upgrade pip' command.
111.png
pysam只支持linux系统下运行，在WIN下缺少DLL文件启动。没办法，只好安装虚拟机和Linux（我用的是ubuntu）在虚拟环境下工作
在准备开发环境前，先试试能不能安装pysam
打开terminal
pip install pysam
Installing collected packages: pysam
Successfully installed pysam-0.15.1
安装成功，OK，终于开始紧锣密鼓的部署开发环境。这里不得不提又来一个坑ubuntu自带pip和python，在我安装好pycharm之后，pycharm无法识别我已经安装好的pysam，我们必须用pip3 install pysam不然他会在原生的python2里面安装。至此，部署环境结束。
任务需求
对BAM文件进行分析，提取所有目标序列，并计算仪器检测上下限GC含量的最大值，最少值。
程序
我直接上最后的程序吧
import pysam
import re
def GC_squence(path_in):
samfile = pysam.AlignmentFile(path_in, "rb")
number = 0
total_GC_MAX = 0 #GC含量最大值計數
total_GC_MIN = 1 #GC含量最小值計數
number_lose= [] #提取失敗序列統計
for line in samfile:
number += 1
if number < 100:
squence = (re.findall("[AGCT][AGCT][AGCT].*[AGCT][AGCT][AGCT]",str(line))) #提取目標序列
print(squence)
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)
else:
number_lose.append(number)
if total_GC > total_GC_MAX:
total_GC_MAX = total_GC
number_MAX =number
if total_GC_MIN > total_GC:
total_GC_MIN = total_GC
number_MIN = number
print("統計總條數%s,最大GC含量%s,編號%s,最小GC含量%s,編號%s,"%(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文件地址
GC_squence(path_in)
试了一下pysam的内置函数，所有内置函数得到结果均为空，但能正常读取序列，最后决定放弃该功能，直接使用正则提取序列。
但是正则提取序列偶尔会出现错误的读长，研究一天依然原因未明。因损失数据为极少量，最后采取跳过策略。
结果
99999 0.6652173913043479
統計總條數1853129,最大GC含量0.8928571428571429,編號85128,最小GC含量0.25,編號63095,
提取失敗序列合計108條,編號爲[1897, 2496, 4247, 4357, 4358, 6860, 11949, 12256, 14237, 15802, 16103, 16112, 16114, 17869, 20134, 23055, 23075, 23102, 23127, 23130, 23135, 23145, 23147, 23152, 23153, 23158, 23161, 23162, 23168, 23169, 23173, 23177, 23179, 23180, 23193, 23198, 23206, 23207, 23210, 23220, 23223, 23231, 23232, 23234, 23239, 23244, 23245, 23247, 23248, 23249, 24120, 30210, 30897, 32018, 36447, 39740, 41929, 45007, 45860, 48468, 48502, 48506, 48518, 48525, 48526, 48579, 48582, 49108, 49330, 52864, 53812, 55594, 55626, 55634, 55653, 55664, 55674, 62145, 62949, 62975, 64321, 64856, 64994, 68642, 69522, 71159, 72266, 73190, 75649, 76752, 80220, 82459, 88235, 88255, 88323, 88638, 90107, 90302, 91539, 91576, 91592, 91594, 95568, 95570, 97141, 98528, 99253, 99274]：
