写在前面
手上小朋友课题做的基因结构注释矫正，为此也喊他练练手，做一些有用的东西出来，今天介绍其中一个。目前，有不少朋友开始重视基因结构注释质量，这是一个好事。我们时长听到关注湿实验的朋友被质量不高的基因结构注释折腾得很惨。我们也很清楚，要做好准确的基因结构注释矫正，往往需要参考各式各样的证据支持，其中，尤其 RNAseq 即转录组数据。但在参考转录组数据时，往往会遇到一个麻烦事：RNAseq 数据太多，而电脑屏幕太小，无法一次看全不同的转录组数据。类似问题，在广西南宁时，有一个师妹提及，至于解决办法，我想了两个：
- 合并所有 BAM 为 1 个BAM，当然这个需要合理设计才能确保合并的BAM足够全面同时又准确；
- 轮转BAM panel，如此就手动多一步操作，可以搞定
对于 2. 轮转 Alignment Tracks，我之前想了一个办法，已经实现过。当然我不太清楚是否有人用上了。
。。。Emmm，忽视外文语法问题，可以看到三个 Tracks 摁一次 Alt+Q 就往上挪一个....如此下来，每个都有机会被看到，如此信息会很全面。
对于 1. 合并 BAM，这块我早前有一个大体想法，简化了想法并给小朋友讲了下，练练手。现在一个月过去了，大体拿到一个我觉得应该能用的脚本。喊他写了一个推文过来，具体如下。
一个合并bam文件的脚本
两个脚本，一个是merge_bam.py，用于合并bam文件,一个是get_state.py，用于评估bam文件。脚本可以通过GitHub下载：genome_annotation_scripts
关于merge_bam.py
为什么写这个脚本，很简单。我们做基因结构注释矫正，往往需要不同组织、不同发育时期的转录组数据，这个时候，会得到许多bam文件。然而，我们不可能把所有的bam文件都导入到GSAman里面辅助基因结构注释矫正，因为电脑屏幕只有这么大。怎么办呢，合并bam文件。这个脚本的原理，就是按照基因所在的位置，划分region，非基因区域以固定长度划分region，再去找到每个region覆盖深度最好的bam文件，将数据写入到merged.bam。
怎么使用merge_bam.py
目前这个脚本只能在linux下面使用，因为pysam在wiodows下面用不了，所以第一步就是要下载pysam
首先创建一个环境（这里我指定了一个python版本，因为我发现在有些版本安装不成功）
conda create -n env_name python=3.10.12
然后用conda下载
conda insatall pysam
## 或者
conda insatall pysam -c bioconda
然后运行python merge_bam.py -h 看看有什么参数
有两个必须传入的参数是输入 input_folder_or_files和 gff_file，分别是输入需要合并的bam文件和一个gff文件，而且需要是跟基因组同一版本的注释文件，不然染色体名称和长度不同，可能会报错。
输入的bam文件可以是放一个文件夹里面(这个脚本会识别这个文件夹里面所有.bam为后缀的文件)，也可以是连续输入n个的bam文件，也可以是文件夹加bam文件，想怎么输入就怎么输入，反正gff文件前面只要是有bam文件或者包含bam文件的文件夹都可以，gff文件放在bam文件后面，位置不可转化，所以最简单的一条命令就是
python merge_bam.py bamfile_folder annotation.gff3
然后就得到一个名为output的文件夹，没有指定输出文件夹默认会创建一个output文件夹，里面有三个文件,merged.bam,merged_sorted.bam和merged_sorted.bam.bai,导入merged_sorted.bam到igv里面就可以开心拉基因了。
可以看到，最上面的那条channel是merged_sorted.bam，从MG、tip、RC中它选择了最好的一个bam文件的reads写入到合并的bam文件里面
关于每个参数是什么意思
-o 是指定一个输出文件夹的名字，默认是创建output文件夹
-r 后面接数字，是指定切分区域的长度，是对于非基因区域而言，因为基因区域不管多长都划为一个区域，默认划分区域长度是5000bp，根据自己基因组大小设置
-s 加上这个参数就是代表对需要合并的bam文件排序和建索引，如果你要合并的bam文件没有索引文件，需要加上这个参数，它会将所有需要合并的bam文件排序和建索引，便且输出到文件夹里面
-t 是使用的cup个数,默认是1个，因为我支持了多进程并行运算，所以-t后面的参数越大，运行起来也会越快
所以使用起来就是
python merge_bam.py bamfile_folder annotation.gff3 -o output_dir -t 4 -r 6000 -s
代表对bamfile_folder文件夹里面所有的文件排序建索引后，以6000bp的长度切分非基因区域，以4个cpu并行运算，结果输出到output_dir里面
关于get_state.py
合并好bam文件后，需要对合并的效果进行评估，这个脚本的功能就是评估bam文件的覆盖率和深度，以及有表达的基因数量
怎么使用get_state.py
看看help信息
这个脚本同样需要输入需要评估的bam文件和gff文件，输入的方法的merge_bam.py一样的，gff文件前面放bam文件或者包含bam文件的文件夹
最简单的命令就是
python get_state.py bamfile_folder annotation.gff3
会输出一个文件output_file，里面包含了每个bam文件的覆盖率，覆盖深度，以及有表达的基因数量
打开output_file可以看到每个bam文件的覆盖率和深度，以及有表达的基因数目，我总共合并了9个bam文件，合并后，merged_sorted.bam对全基因组的覆盖率和有表达的基因数量都明显增加了。
那每个参数的意思是什么呢
-o 代表要输出的文件名字
-t 是代表使用的cpu个数
-c 代表覆盖率，默认是70
-d 代表深度，默认是10
如果你认为一个基因的覆盖率达到70%且覆盖深度有10x说明基因表达了，那所有基因是否表达就可以判断了，这个根据自己的需要来定
那么使用起来就是
python get_state.py bamfile_folder annotation.gff3 -o output.txt -t 4 -c 70 -d 10
意思是对bamfile_folder里面的bam文件统计信息，信息输出到output.txt里面，运行用到的cpu个数是4个，判断一个基因是否表达的标准是覆盖率达到70%且深度有10。
