将二代测序短 reads 比对到参考基因组上,并进行排序和去重复
的等处理,用于后续的变异检测。
part 1 软件和数据准备
软件:bwa samtools picard
文件准备:
参考基因组:genome.fasta
测序数据文件: S1_1.fq.gz 和 S1_2.fq.gz
part 2 构建基因组index
重测序很多软件不能直接读取基因组,需要对基因组构建index才能使用。
```{r setup, include=FALSE}
## samtools 建index
samtools faidx genome.fasta
## bwa 构建 index
bwa index genome.fasta
## picard 构建index
java -jar /opt/picard.jar CreateSequenceDictionary R=genome.fasta
```
#picard 生成:
genome.dict
#samtools 生成:
genome.fasta.fai
#bwa 生成
genome.fasta.amb
genome.fasta.ann
genome.fasta.bwt
genome.fasta.pac
genome.fasta.sa
part 3 比对
```
## step1,read比对并排序
bwa mem -t 2 \ # -t 线程数
-R '@RG\tID:S1\tSM:S1\tPL:illumina' \ #添加read group信息
genome.fasta \ # 参考基因组index路径
S1_1.fq.gz \ # fq1数据路径
S1_2.fq.gz \ # fq2数据路径
2>S1.bwa.log | \ # 比对日志
samtools \
sort -@ 2 \ # 线程数
-m 1G \ # 每个线程内存大小
-o S1.sort.bam - # -o 指定输出文件名称
```
```
## step2, 去除重复
java -Xmx4g -XX:ParallelGCThreads=2 -jar /opt/picard.jar \
MarkDuplicates I=S1.sort.bam \ # 输入文件名称
O=S1.sort.markdup.bam \ # 输出文件名称
CREATE_INDEX=true \ # 创建bam index文件
REMOVE_DUPLICATES=true \ # dup reads是否删除
M=S1.marked_dup_metrics.txt
```
```
##step3, 比对率统计
samtools flagstat S1.sort.bam > S1.sort.bam.flagstat
```
```
##step4, 覆盖度统计
samtools coverage S1.sort.bam > S1.sort.bam.coverage