（7）数据比对过滤

目录

一、为参考基因组建立索引

二、使用BWA-MEM 算法进行mappingzhe

三、用samtools进行下面系列操作

(1)将sam格式转为bam格式

（2）查看bam文件的头部

(3)统计比对结果

（4）提取结果，比对上or没比对上

(5)文件排序

(6)将bam文件转为fastq文件

##bwa mapping到参考基因组

一、为参考基因组建立索引

$ bwa index ref.fa  
#为参考基因组建立索引 (可能稍微有点慢，看参考基因组的大小)
    -a BWT构建算法：bwtsw, is of rb2 [default]，bwtsw适用于较长基因组，另外两个使用于短基因组;
    -p 索引的前缀[same as fasta name]；
    -b bwtsw算法模块长度，与-a bwtsw一起使用，[default 10000000];
$ bwa index Umbilicaria muehlenbergii_genomic.fna

二、使用BWA-MEM 算法进行mappingzhe

$ bwa mem [options] ref.fa reads.fq [mates.fq] 
#使用BWA-MEM 算法进行mapping
##Options:
　　-t INT 线程数，默认是1，增加线程数，会减少运行时间。
　　-M 将 shorter split hits 标记为次优，以兼容 Picard’s markDuplicates 软件。
　　-p 若无此参数：输入文件只有1个，则进行单端比对；若输入文件有2个，则作为paired reads进行比对。若加入此参数：则仅以第1个文件作为输入(会忽略第二个输入序列文件，把第一个文件当做单端测序的数据进行比对)，该文件必须是read1.fq和read2.fa进行reads交叉的数据。
　　-R STR 完整的read group的头部，可以用 '\t' 作为分隔符， 在输出的SAM文件中被解释为制表符TAB. read group 的ID，会被添加到输出文件的每一个read的头部。
　　-T INT 当比对的分值比 INT 小时，不输出该比对结果，这个参数只影响输出的结果，不影响比对的过程。
　　-a 将所有的比对结果都输出，包括 single-end 和 unpaired paired-end的 reads，但是这些比对的结果会被标记为次优。
　　-Y 对数据进行soft clipping, 当错配或者gap数过多比对不上时，会对序列进行切除，这里的切除并只是在比对时去掉这部分序列，最终输出结果中序列还是存在的，所以称为soft

$ bwa mem -t 24 -M ~/Qxy/knowngenome/Umbilicaria_muehlenbergii/Umbilicaria_muehlenbergii_genomic.fna BJ_1.fq BJ_2.fq > BJ_UM_mem.sam 
#一般慢，不能挂nohup运行，后面转格式会报错，所以就慢慢来吧

这个代码意思是将测序的序列和已经建立好索引的Um基因组进行比对

三、用samtools进行下面系列操作

(1)将sam格式转为bam格式

$ samtools view -bS BJ_UM_mem.sam -o BJ_UM_mem.bam -@ 10 
#将sam文件转为bam文件 -b 表示输出bam格式 -S 表示输入sam格式
#超级慢 是不是要加个线程

（2）查看bam文件的头部

$ samtools view -h test.bam | head 
#查看bam文件的头部  SQ：SN表示参考序列的名称，LN表示参考序列的长度 PG：比对时使用的工具指令 HD：VN表示版本，SO表示排序方式
$ samtools view -h BJ_UM_mem.bam | head

(3)统计比对结果

$ samtools flagstat [options] <in.bam> 
#用于给出BAM文件的比对结果 统计reads比对情况
  -@ INT #设置读取文件时要使用的额外线程数
  -O FORMAT # 设置输出格式。FORMAT可以设置为'default'， 'json'或'tsv'来选择默认的，json或标签分隔值输出格式。如果不使用此选项，将选择默认格式
$ samtools flagstat -@ 10 BJ_UM_mem.bam > BJ_UM_mem_flagstat_out.txt 
#并将该结果写入到文本中

（4）提取结果，比对上or没比对上

$ samtools view -b -F 4 BJsample.bam > BJsampleXX.bam 
#提取比对到参考基因组上的数据 -F 表示过滤掉符合FLAG的结果  -4 表示该序列没有比对到参考基因组上

$ samtools view -b -f 4 BJ_UM_mem.bam > BJmem_guolv.bam 
#提取没有比对到参考基因组上的数据 -f 表示只输出FLAG结果 

(5)文件排序

$ samtools sort -@ 10 -o test.bam test.sam 
#对文件进行排序 -o 输出文件名；-@ 后面跟线程数 (稍微看了一下前几行，排不排序没有差别啊)
$ samtools sort -@ 10 -o BJ_UM_map_sort.bam BJ_UM_map.bam

$ samtools view -h BJ_UM_map_sort.bam | head

(6)将bam文件转为fastq文件

$ samtools fastq test.bam > test.fastq 
#将bam文件转化为fastq文件 
$ samtools fastq -@ 10 BJ_UM_map_sort.bam > BJ_UM_map_sort.fastq

后续批量操作也可以写个代码把这些步骤串联起来~