WGS学习全流程笔记
2018-08-03 本文已影响6人
天秤座的机器狗
Part1 数据下载
先去Korean Personal Genome Project下载了编号为KPGP-00001的数据。
先说一下KPGP吧，中文名叫韩国个人基因组计划，这里面的数据都是可以免费下载的。
image
image
nohup wget -c -r -nd -np -k -L -p ftp://ftp.kobic.re.kr/pub/KPGP/2015_release_candidate/WGS/KPGP-00001 1>/dev/null 2>&1 &
简单的说一下这个命令中的几个知识点：
1.nohup
nohup命令可以在你退出帐户之后继续运行相应的进程。nohup就是不挂起的意思( no hang up)
2.&
当在前台运行某个作业时，终端被该作业占据；可以在命令后面加上& 实现后台运行
一般情况下，两者同时使用：
nohup command &
ps all 可以查所有后台任务
这样就不用管它了
下载完成后的数据量如下
image
一般大型的文件的下载需要将其md5文件一起下载，来检验下载的文件是否下载完全
md5sum KPGP-00001_L1_R1.fq.gz >md5tmp1.txt
cat md5tmp1.tx
cat KPGP-00001_L1_R1.fq.gz.md5
如果结果一致就可以判断下载完全
此时的目录应该如下：
image
Part2 测序reads的质控
cat fastqc.sh
#!/bin/bash
for i in $(seq 1 6)
do
fastqc /home/ubuntu/WGS_new/input/KPGP-00001_L${i}_R1.fq.gz /home/ubuntu/WGS_new/input/KPGP-00001_L${i}_R2.fq.gz -o /home/ubuntu/WGS_new/output/fastqc_result && echo "KPGP-00001_L${i} fasqtc is done" >fastqc${i}.log
done
nohup bash fastqc.sh &
由于测序数据有6个lane，而且又是pair-end测序，所以用到multiqc来归纳汇总。
cd ~/WGS_new/output/fastqc_result
multiqc *.zip
这一步之后产生的文件如下
├── input
│ ├── KPGP-00001_L1_R1.fq.gz
│ ├── KPGP-00001_L1_R2.fq.gz
│ ├── KPGP-00001_L2_R1.fq.gz
│ ├── KPGP-00001_L2_R2.fq.gz
│ ├── KPGP-00001_L3_R1.fq.gz
│ ├── KPGP-00001_L3_R2.fq.gz
│ ├── KPGP-00001_L4_R1.fq.gz
│ ├── KPGP-00001_L4_R2.fq.gz
│ ├── KPGP-00001_L5_R1.fq.gz
│ ├── KPGP-00001_L5_R2.fq.gz
│ ├── KPGP-00001_L6_R1.fq.gz
│ ├── KPGP-00001_L6_R2.fq.gz
│ └── nohup.out
├── output
│ └── fastqc_result
│ ├── KPGP-00001_L1_R1_fastqc.html
│ ├── KPGP-00001_L1_R1_fastqc.zip
│ ├── KPGP-00001_L2_R1_fastqc.html
│ ├── KPGP-00001_L2_R1_fastqc.zip
│ ├── KPGP-00001_L2_R2_fastqc.html
│ ├── KPGP-00001_L2_R2_fastqc.zip
│ ├── KPGP-00001_L3_R1_fastqc.html
│ ├── KPGP-00001_L3_R1_fastqc.zip
│ ├── KPGP-00001_L3_R2_fastqc.html
│ ├── KPGP-00001_L3_R2_fastqc.zip
│ ├── KPGP-00001_L4_R1_fastqc.html
│ ├── KPGP-00001_L4_R1_fastqc.zip
│ ├── KPGP-00001_L4_R2_fastqc.html
│ ├── KPGP-00001_L4_R2_fastqc.zip
│ ├── KPGP-00001_L5_R1_fastqc.html
│ ├── KPGP-00001_L5_R1_fastqc.zip
│ ├── KPGP-00001_L5_R2_fastqc.html
│ ├── KPGP-00001_L5_R2_fastqc.zip
│ ├── KPGP-00001_L6_R1_fastqc.html
│ ├── KPGP-00001_L6_R1_fastqc.zip
│ ├── KPGP-00001_L6_R2_fastqc.html
│ ├── KPGP-00001_L6_R2_fastqc.zip
│ ├── multiqc_data
│ │ ├── multiqc_data.json
│ │ ├── multiqc_fastqc.txt
│ │ ├── multiqc_general_stats.txt
│ │ ├── multiqc.log
│ │ └── multiqc_sources.txt
│ └── multiqc_report.html
└── script
├── fastqc1.log
├── fastqc2.log
├── fastqc3.log
├── fastqc4.log
├── fastqc5.log
├── fastqc6.log
├── fastqc.sh
└── nohup.out
此处的质量控制指的是原始数据的质量控制，接下来了解一下质控相关的知识
（此处有借鉴微信公众号：解螺旋的矿工的文章）
说到质量控制，先了解什么是质量值
碱基质量值就是能够用来定量描述碱基好坏程度的一个数值。它该如何才能恰当地描述这个结果呢？我们试想一下，如果测序测得越准确，这个碱基的质量就应该越高；反之，测得越不准确，质量值就应该越低。也就是说可以利用碱基被测错的概率来描述它的质量值，错误率越低，质量值就越高。
假定碱基的测序错误率为p_error，质量值为Q
有公式：
Q = -10log(p_error)
如果该碱基的测序错误率是0.01，那么质量值就是20（俗称Q20），如果是0.001，那么质量值就是30（俗称Q30）。Q20和Q30的比例常常被我们用来评价某次测序结果的好坏，比例越高就越好。
现在一般都是使用Phred33这个体系。
基本上可以从下面几个方面来了解原始数据的质量值：
1.reads各位置的碱基质量分布
在做read质量值分析的时候，FastQC并不单独查看具体某一条read中碱基的质量值，而是将Fastq文件中所有的read数据都综合起来一起分析。如下图，来自KPGP-00001_L1_R1
image.png
这个图的横轴是read上碱基的位置，纵轴是碱基质量值。
我们可以看到read上的每一个位置都有一个黄色的箱型图表示在该位置上所有碱基的质量分布情况。
2.碱基总体质量值分布
fastqc还会提供一个碱基总体的质量分布如下：
image.png
只要大部分都高于20，那么就比较正常
另外，Q20和Q30的比例是我们衡量测序质量的一个重要指标。一般来说，对于二代测序，最好是达到Q20的碱基要在95%以上（最差不低于90%），Q30要求大于85%（最差也不要低于80%）
3.reads各个位置上碱基比例分布
image.png
这个是为了分析碱基的分离程度。何为碱基分离？我们知道AT配对，CG配对，假如测序过程是比较随机的话（随机意味着好），那么在每个位置上A和T比例应该差不多，C和G的比例也应该差不多，如上图所示，两者之间即使有偏差也不应该太大，最好平均在1%以内，如果过高，除非有合理的原因，比如某些特定的捕获测序所致。
4.GC含量分布图
image.png
GC含量指的是G和C这两中碱基占总碱基的比例。二代测序平台或多或少都存在一定的测序偏向性，我们可以通过查看这个值来协助判断测序过程是否足够随机。对于人类来说，我们基因组的GC含量一般在40%左右。因此，如果发现GC含量的图谱明显偏离这个值那么说明测序过程存在较高的序列偏向性，结果就是基因组中某些特定区域被反复测序的几率高于平均水平，除了覆盖度会有偏离之后，将会影响下游的变异检测和CNV分析。图中的GC含量约在39%，处于正常范围。
5.N含量分布图
image.png
N在测序数据中一般是不应该出现的，如果出现则意味着，测序的光学信号无法被清晰分辨，如果这种情况多的话，往往意味着测序系统或者测序试剂的错误。图中接近0%。
6.接头序列
image.png
在测序之前需要构建测序文库，测序接头就是在这个时候加上的，其目的一方面是为了能够结合到flowcell上，另一方面是当有多个样本同时测序的时候能够利用接头信息进行区分。当测序read的长度大于被测序的DNA片段时，就会在read的末尾测到这些接头序列。一般的WGS测序是不会测到这些接头序列的，因为构建WGS测序的文库序列（插入片段）都比较长，约几百bp，而read的测序长度都在100bp-150bp这个范围。不过在进行一些RNA测序的时候，由于它们的序列本来就比较短很多只有几十bp长（特别是miRNA），那么就很容易会出现read测通的现象，这个时候就会在read的末尾测到这些接头序列。
image.png
最后，如果测序质量很差的情况下，则需要切除测序接头序列和read的低质量序列
工具有:SOAPnuke、cutadapt、untrimmed等,但是，常用的一般是SOAPnuke、cutadapt、untrimmed等。
Part3 比对（mapping）
step1
mkdir /home/ubuntu/WGS_new/input/genome
cd /home/ubuntu/WGS_new/input/genome
wget http://hgdownload.cse.ucsc.edu/goldenPath/hg38/bigZips/hg38.fa.gz
#对下载的基因组建立索引，下面是建索引的脚本
#!/bin/bash
bwa index -a bwtsw -p /home/ubuntu/WGS_new/input/genome/hg38 /home/ubuntu/WGS_new/input/genome/Homo_sapiens_assembly38.fasta && echo "bwa index of hg38 is done"
#运行脚本
nohup bash bwa_index.sh &
索引构建完毕后的目录如下：
├── input
│ ├── genome
│ │ ├── hg38.amb
│ │ ├── hg38.ann
│ │ ├── hg38.bwt
│ │ ├── hg38.fa
│ │ ├── hg38.pac
│ │ ├── hg38.sa
│ │ └── nohup.out
│ ├── KPGP-00001_L1_R1.fq.gz
│ ├── KPGP-00001_L1_R2.fq.gz
│ ├── KPGP-00001_L2_R1.fq.gz
│ ├── KPGP-00001_L2_R2.fq.gz
│ ├── KPGP-00001_L3_R1.fq.gz
│ ├── KPGP-00001_L3_R2.fq.gz
│ ├── KPGP-00001_L4_R1.fq.gz
│ ├── KPGP-00001_L4_R2.fq.gz
│ ├── KPGP-00001_L5_R1.fq.gz
│ ├── KPGP-00001_L5_R2.fq.gz
│ ├── KPGP-00001_L6_R1.fq.gz
│ ├── KPGP-00001_L6_R2.fq.gz
│ └── nohup.out
├── output
│ └── fastqc_result
│ ├── KPGP-00001_L1_R1_fastqc.html
│ ├── KPGP-00001_L1_R1_fastqc.zip
│ ├── KPGP-00001_L2_R1_fastqc.html
│ ├── KPGP-00001_L2_R1_fastqc.zip
│ ├── KPGP-00001_L2_R2_fastqc.html
│ ├── KPGP-00001_L2_R2_fastqc.zip
│ ├── KPGP-00001_L3_R1_fastqc.html
│ ├── KPGP-00001_L3_R1_fastqc.zip
│ ├── KPGP-00001_L3_R2_fastqc.html
│ ├── KPGP-00001_L3_R2_fastqc.zip
│ ├── KPGP-00001_L4_R1_fastqc.html
│ ├── KPGP-00001_L4_R1_fastqc.zip
│ ├── KPGP-00001_L4_R2_fastqc.html
│ ├── KPGP-00001_L4_R2_fastqc.zip
│ ├── KPGP-00001_L5_R1_fastqc.html
│ ├── KPGP-00001_L5_R1_fastqc.zip
│ ├── KPGP-00001_L5_R2_fastqc.html
│ ├── KPGP-00001_L5_R2_fastqc.zip
│ ├── KPGP-00001_L6_R1_fastqc.html
│ ├── KPGP-00001_L6_R1_fastqc.zip
│ ├── KPGP-00001_L6_R2_fastqc.html
│ ├── KPGP-00001_L6_R2_fastqc.zip
│ ├── multiqc_data
│ │ ├── multiqc_data.json
│ │ ├── multiqc_fastqc.txt
│ │ ├── multiqc_general_stats.txt
│ │ ├── multiqc.log
│ │ └── multiqc_sources.txt
│ └── multiqc_report.html
└── script
├── bwa_index.sh
├── fastqc1.log
├── fastqc2.log
├── fastqc3.log
├── fastqc4.log
├── fastqc5.log
├── fastqc6.log
├── fastqc.sh
└── nohup.out
step2
cd /home/ubuntu/WGS_new/output
mkdir /home/ubuntu/WGS_new/output/bwa_mapping
cd /home/ubuntu/WGS_new/script
#比对的脚本如下
#!/bin/bash
bwa mem -t 4 -R '@RG\tID:KPGP-00001_L1\tPL:illumina\tSM:KPGP-00001' /home/ubuntu/WGS_new/input/genome/hg38 /home/ubuntu/WGS_new/input/KPGP-00001_L1_R1.fq.gz /home/ubuntu/WGS_new/input/KPGP-00001_L1_R2.fq.gz >/home/ubuntu/WGS_new/output/bwa_mapping/KPGP-00001_L1.sam && echo "KPGP-00001_L1 mapping is done"
bwa mem -t 4 -R '@RG\tID:KPGP-00001_L2\tPL:illumina\tSM:KPGP-00001' /home/ubuntu/WGS_new/input/genome/hg38 /home/ubuntu/WGS_new/input/KPGP-00001_L2_R1.fq.gz /home/ubuntu/WGS_new/input/KPGP-00001_L2_R2.fq.gz >/home/ubuntu/WGS_new/output/bwa_mapping/KPGP-00001_L2.sam && echo "KPGP-00001_L2 mapping is done"
bwa mem -t 4 -R '@RG\tID:KPGP-00001_L3\tPL:illumina\tSM:KPGP-00001' /home/ubuntu/WGS_new/input/genome/hg38 /home/ubuntu/WGS_new/input/KPGP-00001_L3_R1.fq.gz /home/ubuntu/WGS_new/input/KPGP-00001_L3_R2.fq.gz >/home/ubuntu/WGS_new/output/bwa_mapping/KPGP-00001_L3.sam && echo "KPGP-00001_L3 mapping is done"
bwa mem -t 4 -R '@RG\tID:KPGP-00001_L4\tPL:illumina\tSM:KPGP-00001' /home/ubuntu/WGS_new/input/genome/hg38 /home/ubuntu/WGS_new/input/KPGP-00001_L4_R1.fq.gz /home/ubuntu/WGS_new/input/KPGP-00001_L4_R2.fq.gz >/home/ubuntu/WGS_new/output/bwa_mapping/KPGP-00001_L4.sam && echo "KPGP-00001_L4 mapping is done"
bwa mem -t 4 -R '@RG\tID:KPGP-00001_L5\tPL:illumina\tSM:KPGP-00001' /home/ubuntu/WGS_new/input/genome/hg38 /home/ubuntu/WGS_new/input/KPGP-00001_L5_R1.fq.gz /home/ubuntu/WGS_new/input/KPGP-00001_L5_R2.fq.gz >/home/ubuntu/WGS_new/output/bwa_mapping/KPGP-00001_L5.sam && echo "KPGP-00001_L5 mapping is done"
bwa mem -t 4 -R '@RG\tID:KPGP-00001_L6\tPL:illumina\tSM:KPGP-00001' /home/ubuntu/WGS_new/input/genome/hg38 /home/ubuntu/WGS_new/input/KPGP-00001_L6_R1.fq.gz /home/ubuntu/WGS_new/input/KPGP-00001_L6_R2.fq.gz >/home/ubuntu/WGS_new/output/bwa_mapping/KPGP-00001_L6.sam && echo "KPGP-00001_L6 mapping is done"
#运行脚本
nohup bash bwa_mapping.sh &
比对后生成sam文件，将sam文件转换成二进制格式的bam文件
#转换的脚本如下
#!/bin/bash
for i in $(seq 1 6)
do
samtools view -Sb /home/ubuntu/WGS_new/output/bwa_mapping/KPGP-00001_L${i}.sam >/home/ubuntu/WGS_new/output/bwa_mapping/KPGP-00001_L${i}.bam && echo "KPGP-00001_L${i} samt2bam is done"
done
#运行脚本
nohup bash sam2bam.sh &
#运行结束后，再把sam格式文件删除，其实可以不删除，但是太占储存空间了
cd /home/ubuntu/WGS_new/output/bwa_mapping
rm *.sam
结果如下：
├── input
│ ├── genome
│ │ ├── hg38.amb
│ │ ├── hg38.ann
│ │ ├── hg38.bwt
│ │ ├── hg38.fa
│ │ ├── hg38.pac
│ │ ├── hg38.sa
│ │ └── nohup.out
│ ├── KPGP-00001_L1_R1.fq.gz
│ ├── KPGP-00001_L1_R2.fq.gz
│ ├── KPGP-00001_L2_R1.fq.gz
│ ├── KPGP-00001_L2_R2.fq.gz
│ ├── KPGP-00001_L3_R1.fq.gz
│ ├── KPGP-00001_L3_R2.fq.gz
│ ├── KPGP-00001_L4_R1.fq.gz
│ ├── KPGP-00001_L4_R2.fq.gz
│ ├── KPGP-00001_L5_R1.fq.gz
│ ├── KPGP-00001_L5_R2.fq.gz
│ ├── KPGP-00001_L6_R1.fq.gz
│ ├── KPGP-00001_L6_R2.fq.gz
│ └── nohup.out
├── output
│ ├── bwa_mapping
│ │ ├── KPGP-00001_L1.bam
│ │ ├── KPGP-00001_L2.bam
│ │ ├── KPGP-00001_L3.bam
│ │ ├── KPGP-00001_L4.bam
│ │ ├── KPGP-00001_L5.bam
│ │ └── KPGP-00001_L6.bam
│ └── fastqc_result
│ ├── KPGP-00001_L1_R1_fastqc.html
│ ├── KPGP-00001_L1_R1_fastqc.zip
│ ├── KPGP-00001_L2_R1_fastqc.html
│ ├── KPGP-00001_L2_R1_fastqc.zip
│ ├── KPGP-00001_L2_R2_fastqc.html
│ ├── KPGP-00001_L2_R2_fastqc.zip
│ ├── KPGP-00001_L3_R1_fastqc.html
│ ├── KPGP-00001_L3_R1_fastqc.zip
│ ├── KPGP-00001_L3_R2_fastqc.html
│ ├── KPGP-00001_L3_R2_fastqc.zip
│ ├── KPGP-00001_L4_R1_fastqc.html
│ ├── KPGP-00001_L4_R1_fastqc.zip
│ ├── KPGP-00001_L4_R2_fastqc.html
│ ├── KPGP-00001_L4_R2_fastqc.zip
│ ├── KPGP-00001_L5_R1_fastqc.html
│ ├── KPGP-00001_L5_R1_fastqc.zip
│ ├── KPGP-00001_L5_R2_fastqc.html
│ ├── KPGP-00001_L5_R2_fastqc.zip
│ ├── KPGP-00001_L6_R1_fastqc.html
│ ├── KPGP-00001_L6_R1_fastqc.zip
│ ├── KPGP-00001_L6_R2_fastqc.html
│ ├── KPGP-00001_L6_R2_fastqc.zip
│ ├── multiqc_data
│ │ ├── multiqc_data.json
│ │ ├── multiqc_fastqc.txt
│ │ ├── multiqc_general_stats.txt
│ │ ├── multiqc.log
│ │ └── multiqc_sources.txt
│ └── multiqc_report.html
└── script
├── bwa_index.sh
├── bwa_mapping.sh
├── fastqc1.log
├── fastqc2.log
├── fastqc3.log
├── fastqc4.log
├── fastqc5.log
├── fastqc6.log
├── fastqc.sh
├── nohup.out
└── sam2bam.sh
step3
对生成的bam文件按position排序
#脚本如下
#!/bin/bash
for i in $(seq 1 6)
do
samtools sort -@ 4 -m 4G -O bam -o /home/ubuntu/WGS_new/output/bwa_mapping/KPGP-00001_L${i}.sorted.bam /home/ubuntu/WGS_new/output/bwa_mapping/KPGP-00001_L${i}.bam && echo "/KPGP-00001_L${i}.sorted.bam is done"
done
此时所有结果如下
├── input
│ ├── genome
│ │ ├── hg38.amb
│?? │?? ├── hg38.ann
│?? │?? ├── hg38.bwt
│?? │?? ├── hg38.fa
│?? │?? ├── hg38.pac
│?? │?? ├── hg38.sa
│?? │?? └── nohup.out
│?? ├── KPGP-00001_L1_R1.fq.gz
│?? ├── KPGP-00001_L1_R2.fq.gz
│?? ├── KPGP-00001_L2_R1.fq.gz
│?? ├── KPGP-00001_L2_R2.fq.gz
│?? ├── KPGP-00001_L3_R1.fq.gz
│?? ├── KPGP-00001_L3_R2.fq.gz
│?? ├── KPGP-00001_L4_R1.fq.gz
│?? ├── KPGP-00001_L4_R2.fq.gz
│?? ├── KPGP-00001_L5_R1.fq.gz
│?? ├── KPGP-00001_L5_R2.fq.gz
│?? ├── KPGP-00001_L6_R1.fq.gz
│?? ├── KPGP-00001_L6_R2.fq.gz
│?? └── nohup.out
├── output
│?? ├── bwa_mapping
│?? │?? ├── index
│?? │?? ├── KPGP-00001_L1.bam
│?? │?? ├── KPGP-00001_L1.sorted.bam
│?? │?? ├── KPGP-00001_L2.bam
│?? │?? ├── KPGP-00001_L2.sorted.bam
│?? │?? ├── KPGP-00001_L3.bam
│?? │?? ├── KPGP-00001_L3.sorted.bam
│?? │?? ├── KPGP-00001_L4.bam
│?? │?? ├── KPGP-00001_L4.sorted.bam
│?? │?? ├── KPGP-00001_L5.bam
│?? │?? ├── KPGP-00001_L5.sorted.bam
│?? │?? ├── KPGP-00001_L6.bam
│?? │?? └── KPGP-00001_L6.sorted.bam
│?? └── fastqc_result
│?? ├── KPGP-00001_L1_R1_fastqc.html
│?? ├── KPGP-00001_L1_R1_fastqc.zip
│?? ├── KPGP-00001_L2_R1_fastqc.html
│?? ├── KPGP-00001_L2_R1_fastqc.zip
│?? ├── KPGP-00001_L2_R2_fastqc.html
│?? ├── KPGP-00001_L2_R2_fastqc.zip
│?? ├── KPGP-00001_L3_R1_fastqc.html
│?? ├── KPGP-00001_L3_R1_fastqc.zip
│?? ├── KPGP-00001_L3_R2_fastqc.html
│?? ├── KPGP-00001_L3_R2_fastqc.zip
│?? ├── KPGP-00001_L4_R1_fastqc.html
│?? ├── KPGP-00001_L4_R1_fastqc.zip
│?? ├── KPGP-00001_L4_R2_fastqc.html
│?? ├── KPGP-00001_L4_R2_fastqc.zip
│?? ├── KPGP-00001_L5_R1_fastqc.html
│?? ├── KPGP-00001_L5_R1_fastqc.zip
│?? ├── KPGP-00001_L5_R2_fastqc.html
│?? ├── KPGP-00001_L5_R2_fastqc.zip
│?? ├── KPGP-00001_L6_R1_fastqc.html
│?? ├── KPGP-00001_L6_R1_fastqc.zip
│?? ├── KPGP-00001_L6_R2_fastqc.html
│?? ├── KPGP-00001_L6_R2_fastqc.zip
│?? ├── multiqc_data
│?? │?? ├── multiqc_data.json
│?? │?? ├── multiqc_fastqc.txt
│?? │?? ├── multiqc_general_stats.txt
│?? │?? ├── multiqc.log
│?? │?? └── multiqc_sources.txt
│?? └── multiqc_report.html
└── script
├── bam_index.sh
├── bwa_index.sh
├── bwa_mapping.sh
├── fastqc1.log
├── fastqc2.log
├── fastqc3.log
├── fastqc4.log
├── fastqc5.log
├── fastqc6.log
├── fastqc.sh
├── nohup.out
└── sam2bam.sh
以上的步骤就完成了比对的过程
在此基础上，做一些以下的横向探索:
1.分别利用 samtools和 bedtools 对测序深度和覆盖度进行统计，并可视化
2.把bam文件按染色体分成小文件，并用IGV查看
3.bam格式转化为bw格式看测序深度分布
4.测序深度和GC含量的关系
1.分别利用 samtools和 bedtools 对测序深度和覆盖度进行统计，并可视化
用samtools 对测序深度进行统计
#!/bin/bash
for i in $(seq 1 6)
do
samtools depth /home/ubuntu/WGS_new/output/bwa_mapping/KPGP-00001_L${i}.sorted.bam 1>/home/ubuntu/WGS_new/output/bwa_mapping/samtools_depth/KPGP-00001_L${i}.sorted.bam.depth 2>&1 && echo "KPGP-00001_L${i}.sorted.bam depth count is done"
done
这个软件得到的文件格式如下
chr1 9999 4
chr1 10000 5
chr1 10001 51
chr1 10002 161
chr1 10003 250
chr1 10004 278
chr1 10005 365
chr1 10006 390
chr1 10007 391
chr1 10008 397
chr1 10009 406
chr1 10010 407
chr1 10011 414
chr1 10012 415
chr1 10013 417
chr1 10014 420
chr1 10015 420
chr1 10016 419
chr1 10017 421
chr1 10018 422
chr1 10019 425
chr1 10020 427
chr1 10021 430
chr1 10022 432
chr1 10023 435
chr1 10024 435
chr1 10025 437
chr1 10026 439
chr1 10027 439
chr1 10028 440
chr1 10029 440
chr1 10030 442
...
...
第一列是染色体号，第二列是对应的单个碱基位置，第三列是单个碱基的深度
需要进一步统计到每个染色体对应的测序深度
这一步的python脚本如下
import sys
args=sys.argv
filename=args[1]
SumDict={}
for line in open(filename):
lineL=line.strip().split("\t")
ID=lineL[0]
depth=int(lineL[2])
if ID not in SumDict:
SumDict[ID]=depth
else:
SumDict[ID]+=depth
for k,v in SumDict.items():
print(k,v,sep="\t")
-----------------------------
nohup time python3 sumDepth.py /home/ubuntu/WGS_new/output/bwa_mapping/samtools_depth/KPGP-00001_L1.sorted.bam.depth > /home/ubuntu/WGS_new/output/bwa_mapping/samtools_depth/KPGP-00001_L1.SumOfDepth.txt &
nohup time python3 sumDepth.py /home/ubuntu/WGS_new/output/bwa_mapping/samtools_depth/KPGP-00001_L2.sorted.bam.depth > /home/ubuntu/WGS_new/output/bwa_mapping/samtools_depth/KPGP-00001_L2.SumOfDepth.txt &
nohup time python3 sumDepth.py /home/ubuntu/WGS_new/output/bwa_mapping/samtools_depth/KPGP-00001_L3.sorted.bam.depth > /home/ubuntu/WGS_new/output/bwa_mapping/samtools_depth/KPGP-00001_L3.SumOfDepth.txt &
nohup time python3 sumDepth.py /home/ubuntu/WGS_new/output/bwa_mapping/samtools_depth/KPGP-00001_L4.sorted.bam.depth > /home/ubuntu/WGS_new/output/bwa_mapping/samtools_depth/KPGP-00001_L4.SumOfDepth.txt &
nohup time python3 sumDepth.py /home/ubuntu/WGS_new/output/bwa_mapping/samtools_depth/KPGP-00001_L5.sorted.bam.depth > /home/ubuntu/WGS_new/output/bwa_mapping/samtools_depth/KPGP-00001_L5.SumOfDepth.txt &
nohup time python3 sumDepth.py /home/ubuntu/WGS_new/output/bwa_mapping/samtools_depth/KPGP-00001_L6.sorted.bam.depth > /home/ubuntu/WGS_new/output/bwa_mapping/samtools_depth/KPGP-00001_L6.SumOfDepth.txt &
得到的统计结果还需要把带“_”部分除掉
#!/bin/bash
for i in $(seq 1 6)
do
grep -v "_" /home/ubuntu/WGS_new/output/bwa_mapping/samtools_depth/KPGP-00001_L${i}.SumOfDepth.txt >/home/ubuntu/WGS_new/output/bwa_mapping/samtools_depth/KPGP-00001_L${i}.SumOfDepth.final.txt
done
然后将所有的生成文件排序后合并，脚本如下
cd /home/ubuntu/WGS_new/output/bwa_mapping/samtools_depth
#脚本
--------
#!/bin/bash
for i in $(seq 1 6)
do
sed -i "s/chr//" KPGP-00001_L${i}.SumOfDepth.final.txt #首次去掉chr，方便排序
sort -k1,1n KPGP-00001_L${i}.SumOfDepth.final.txt > KPGP-00001_L${i}.SumOfDepth.sort.final.txt #对第一列按数值排序
sed -i 's/^/chr/' KPGP-00001_L${i}.SumOfDepth.sort.final.txt #再将chr添加回来
done
paste KPGP-00001_L*.SumOfDepth.sort.final.txt >KPGP-00001_all.SumOfDepth.sort.final.txt
awk '{sum=$2+$$4+$6+$8+$10+$12;print $1,sum;}' KPGP-00001_all.SumOfDepth.sort.final.txt >KPGP-00001.SumOfDepth.final_count.txt
---------
#运行脚本
nohup bash paste_sort.sh &
-----------
sed -i '1 i Chrom\tsum_of_depth' KPGP-00001.SumOfDepth.final_count.txt
cat KPGP-00001.SumOfDepth.final_count.txt
-----------
Chrom sum_of_depth
chrM 34591025
chrX 4093132439
chrY 125469175
chr1 6482954037
chr2 6376037101
chr3 5286606057
chr4 5330867456
chr5 4949134249
chr6 4420827669
chr7 4173524521
chr8 3842370943
chr9 3228529809
chr10 3626998764
chr11 3548159595
chr12 3490214776
chr13 2618593673
chr14 2359781279
chr15 2098133774
chr16 2341297030
chr17 1911455487
chr18 2082861061
chr19 1341657686
chr20 1741460772
chr21 1203547746
chr22 988334122
sed -i 's/ /\t/' KPGP-00001.SumOfDepth.final_count.txt
以上就将所有的染色体的覆盖位点的测序深度总和统计出来了
接下来计算覆盖度
因为samtools depth计算的是所有覆盖的每个碱基的深度，所以将对应染色体的统计信息提取，并计算有多少个碱基，得到的碱基数就是对应的覆盖长度。
所以有如下脚本统计
#!/bin/bash
for i in $(seq 1 22)
do
grep "chr${i} " /home/ubuntu/WGS_new/output/bwa_mapping/samtools_depth/KPGP-00001_L1.sorted.bam.depth |wc -l >/home/ubuntu/WGS_new/output/bwa_mapping/samtools_depth/chr${i}_cov.txt
sed -i "1 i chr${i}" /home/ubuntu/WGS_new/output/bwa_mapping/samtools_depth/chr${i}_cov.txt
done
grep "chrX " /home/ubuntu/WGS_new/output/bwa_mapping/samtools_depth/KPGP-00001_L1.sorted.bam.depth |wc -l >/home/ubuntu/WGS_new/output/bwa_mapping/samtools_depth/chrX_cov.txt #x染色体
sed -i "1 i chrX" /home/ubuntu/WGS_new/output/bwa_mapping/samtools_depth/chrX_cov.txt
grep "chrY " /home/ubuntu/WGS_new/output/bwa_mapping/samtools_depth/KPGP-00001_L1.sorted.bam.depth |wc -l >/home/ubuntu/WGS_new/output/bwa_mapping/samtools_depth/chrY_cov.txt #y染色体
sed -i "1 i chrY" /home/ubuntu/WGS_new/output/bwa_mapping/samtools_depth/chrY_cov.txt
grep "chrM " /home/ubuntu/WGS_new/output/bwa_mapping/samtools_depth/KPGP-00001_L1.sorted.bam.depth |wc -l >/home/ubuntu/WGS_new/output/bwa_mapping/samtools_depth/chrM_cov.txt #线粒体
sed -i "1 i chrM" /home/ubuntu/WGS_new/output/bwa_mapping/samtools_depth/chrM_cov.txt
这里计算的是一个lane的coverage，应该要计算所有6个lane的coverage，再求mean值，或者提前将6个lane的bam文件merge后再处理。此处由于系统存储空间问题，只对一个lane作简单统计。
sort -m chr*_cov.txt|sed -n 1~2p> chr_cov1.txt
sort -m chr*_cov.txt|sed -n 2~2p> chr_cov2.txt
paste chr_cov*.txt>chr_cov_all.txt
sed 's/chr//' chr_cov_all.txt|sort -k1,1n|sed 's/^/chr/'|sed '1 i Chrom\tcoverage' >chr_cov.txt
cat chr_cov.txt
Chrom coverage
chrM 16569
chrX 151545600
chrY 4004642
chr1 224300366
chr2 235466712
chr3 194391189
chr4 185704231
chr5 176555626
chr6 164731880
chr7 155128428
chr8 141501464
chr9 118343356
chr10 129732215
chr11 131193733
chr12 130177621
chr13 95749291
chr14 87847379
chr15 81867181
chr16 78578328
chr17 78271424
chr18 77398752
chr19 54075815
chr20 62148856
chr21 38164294
chr22 36926516
计算参考序列染色体的长度
这里首先需要将hg38中的重复的N序列去除
首先想到的是
grep -v "N" hg38.fa >hg38.cut.fa
这里有一个小问题，就是如果N藏在正常ATGC序列之中的话，就很容易去除掉正常的碱基
所以改成sed命令
sed 's/N//g' hg38.fa >hg38.sed.fa #把所有的N替换成空
sed '/^$/d' hg38.sed.fa >hg38.final.fa #去除空行
sed是将N替换成空，并不是将含N的行除去
虽然两个命令得到结果一样，但还是选择后者是正确的。
写了2个python脚本分别统计1.hg38原始文件对应各条染色体的长度
2.hg38去除N之后的各条染色体的长度
3.hg38原始文件的各条染色体的N含量
脚本如下
import sys
args=sys.argv
filename=args[1]
seq_id={}
for line in open(filename):
line=line.strip()
if line.startswith(">"):
seq_name = line
if seq_name not in seq_id:
seq_id[seq_name]=''
else:
seq_id[seq_name]+=line.strip()
for k,v in seq_id.items():
seq_length=len(v)
name=k[1:]
print(name,seq_length,sep="\t")
---------------------------------------------------
#运行脚本
--------------------
nohup python3 count_length.py hg38.fa > hg38.chr_length.txt &
nohup python3 count_length.py hg38.final.fa > hg38.cutN.chr_length.txt &
-----------------------------------
import sys
args=sys.argv
filename=args[1]
seq_id={}
for line in open(filename):
line=line.strip()
if line.startswith(">"):
seq_name = line
if seq_name not in seq_id:
seq_id[seq_name]=''
else:
seq_id[seq_name]+=line
for k,v in seq_id.items():
N_count=v.count("N")
name=k[1:]
print(k,N_count,sep="\t")
--------------------------
#运行脚本
nohup python3 count_N.py hg38.fa >hg38.chr.countN.txt &
这两个脚本本身没有问题，而且对于fasta格式读取，把ID 和序列先构建字典的方式也是之前总结出来的
但是！！！！
对于大型的文件，先构建字典，在读取字典的过程会花费太多的时间，所以这两个脚本运行了如下的时间，程序还没有跑完
image.png
所以，需要修改代码，不要用这种构建字典的方式。
修改后的脚本如下：
#统计长度的脚本
import sys
args=sys.argv
filename=args[1]
seq_length=0
for line in open(filename):
line=line.strip()
if line.startswith(">"):
if seq_length:
print(seq_length)
seq_length=0
print(line)
else:
seq_length+=len(line)
print(seq_length)
--------------------------------------
#统计N数量的脚本
import sys
args=sys.argv
filename=args[1]
num_N=0
for line in open(filename):
line=line.strip()
if line.startswith(">"):
if num_N :
print (num_N)
num_N=0
print (line)
else:
num_N+=line.count("N")
print(num_N)
-----------
优化后的代码只花了2到3分钟就完成任务，可见代码优化的重要性。
因为得到结果不仅仅有chr1~22+x+y+M的信息，还有其他信息，所以需要提取
提取脚本如下
#!/bin/bash
for i in $(seq 1 22)
do
grep -A 1 ">chr${i}$" hg38.cutN.chr_length.txt > hg38.cutN.chr_length.final${i}.txt
grep -A 1 ">chr${i}$" count_N_new.txt > count_N_new.final${i}.txt
done
#!/bin/bash
for i in $(seq 1 22)
do
grep -A 1 ">chr${i}$" hg38.chr_length.txt >hg38.chr_length.final${i}.txt #提取>chr${i}$以及其下面一行，-A是after的意思，之所以是chr${i}+一个$是因为chr${i}后面加空行
才是想要的结果，如果不清楚可以cat -A 一下，再判断
grep -A 1 ">chr${i}$" hg38.cutN.chr_length.txt > hg38.cutN.chr_length.final${i}.txt
grep -A 1 ">chr${i}$" count_N_new.txt > count_N_new.final${i}.txt
done
grep -A 1 ">chrX$" hg38.chr_length.txt >hg38.chr_length.finalX.txt
grep -A 1 ">chrX$" hg38.cutN.chr_length.txt > hg38.cutN.chr_length.finalX.txt
grep -A 1 ">chrX$" count_N_new.txt > count_N_new.finalX.txt
grep -A 1 ">chrY$" hg38.chr_length.txt >hg38.chr_length.finalY.txt
grep -A 1 ">chrY$" hg38.cutN.chr_length.txt > hg38.cutN.chr_length.finalY.txt
grep -A 1 ">chrY$" count_N_new.txt > count_N_new.finalY.txt
grep -A 1 ">chrM$" hg38.chr_length.txt >hg38.chr_length.finalM.txt
grep -A 1 ">chrM$" hg38.cutN.chr_length.txt > hg38.cutN.chr_length.finalM.txt
grep -A 1 ">chrM$" count_N_new.txt > count_N_new.finalM.txt
rm hg38.chr_length.txt
rm hg38.cutN.chr_length.txt
rm count_N_new.txt
sort -m hg38.chr_length.final*.txt >hg38.chr_length.final_all.txt #将所有hg38.chr_length.final*.txt 文件合并，区别于paste，这个合并成列
sort -m hg38.cutN.chr_length.final*.txt >hg38.cutN.chr_length.final_all.txt
sort -m count_N_new.final*.txt >count_N_new.final_all.txt
sed -n 1~2p hg38.chr_length.final_all.txt >hg38.chr_length.final_all_1.txt #提取奇数行，1~2p中1是第一行，2p是从第一行开始，每2行输出，也就是输出奇数行
sed -n 2~2p hg38.chr_length.final_all.txt >hg38.chr_length.final_all_2.txt #提取偶数行，2~2p中1是第一行，2p是从第一行开始，每2行输出，也就是输出偶数行
paste hg38.chr_length.final_all_1.txt hg38.chr_length.final_all_2.txt > hg38.chr_length.final_all_3.txt #hg38每条染色体的长度，包括N
这里包括3个知识点：
- grep -A 1 string file 提取file中匹配到string的行和下面一行，1可以改成其他数字，表示多行
grep -B 1 string file 上面
grep -C 1 string file 上面和下面
2.cat -A 之后再找匹配的字符串比较准确
- sort -m 将文件合并成列，区别与paste
- sed -n 1~2p 提取奇数行
sed -n 2~2p 提取偶数行
当然awk也很容易实现
上面提取之后的文件排序，加列标题，然后合并
脚本如下：
sed 's/>chr//' count_N_new.final_all_3.txt |sort -k1,1n |sed 's/^/chr/'|sed '1 i Chrom\tN_Number' >N_number.txt #排序，然后添加列标题
sed 's/>chr//' hg38.cutN.chr_length.final_all_3.txt |sort -k1,1n |sed 's/^/chr/'|sed '1 i Chrom\tCut_N_length' >hg38.cutN.len.txt
sed 's/>chr//' hg38.chr_length.final_all_3.txt |sort -k1,1n |sed 's/^/chr/'|sed '1 i Chrom\tlength' >hg38.len.txt
paste hg38.len.txt hg38.cutN.len.txt N_number.txt |cut -f 1,2,4,6 >hg38.count_all_final.txt
less -S hg38.count_all_final.txt
-------------------------------------------------
Chrom length Cut_N_length N_Number
chrM 16569 16568 1
chrX 156040895 154893130 1147765
chrY 57227415 26415183 30812232
chr1 248956422 230481193 18475229
chr2 242193529 240548238 1645291
chr3 198295559 198100139 195420
chr4 190214555 189752667 461888
chr5 181538259 181265378 272881
chr6 170805979 170078724 727255
chr7 159345973 158970132 375841
chr8 145138636 144768136 370500
chr9 138394717 121790590 16604127
chr10 133797422 133263135 534287
chr11 135086622 134533742 552880
chr12 133275309 133138039 137270
chr13 114364328 97983128 16381200
chr14 107043718 90568149 16475569
chr15 101991189 84641348 17349841
chr16 90338345 81805944 8532401
chr17 83257441 82921074 336367
chr18 80373285 80089658 283627
chr19 58617616 58440758 176858
chr20 64444167 63944581 499586
chr21 46709983 40088683 6621300
chr22 50818468 39159777 11658691
--------------
结果一目了然
-------------------------------------------------------------------------------
最后再把测序深度、覆盖度、参考序列长度合并
cd /home/ubuntu/WGS_new/output/bwa_mapping/samtools_depth/new
ls -lh
------
total 12K
-rw-rw-r-- 1 ubuntu ubuntu 387 Aug 7 22:11 chr_cov.txt
-rw-rw-r-- 1 ubuntu ubuntu 834 Aug 7 21:51 hg38.count_all_final.txt
-rw-rw-r-- 1 ubuntu ubuntu 429 Aug 7 22:18 KPGP-00001.SumOfDepth.final_count.txt
-------
paste * |cut -f 1,2,4,5,6,8 >count_all.txt
cat count_all.txt
Chrom coverage length Cut_N_length N_Number sum_of_depth
chrM 16569 16569 16568 1 34591025
chrX 151545600 156040895 154893130 1147765 4093132439
chrY 4004642 57227415 26415183 30812232 125469175
chr1 224300366 248956422 230481193 18475229 6482954037
chr2 235466712 242193529 240548238 1645291 6376037101
chr3 194391189 198295559 198100139 195420 5286606057
chr4 185704231 190214555 189752667 461888 5330867456
chr5 176555626 181538259 181265378 272881 4949134249
chr6 164731880 170805979 170078724 727255 4420827669
chr7 155128428 159345973 158970132 375841 4173524521
chr8 141501464 145138636 144768136 370500 3842370943
chr9 118343356 138394717 121790590 16604127 3228529809
chr10 129732215 133797422 133263135 534287 3626998764
chr11 131193733 135086622 134533742 552880 3548159595
chr12 130177621 133275309 133138039 137270 3490214776
chr13 95749291 114364328 97983128 16381200 2618593673
chr14 87847379 107043718 90568149 16475569 2359781279
chr15 81867181 101991189 84641348 17349841 2098133774
chr16 78578328 90338345 81805944 8532401 2341297030
chr17 78271424 83257441 82921074 336367 1911455487
chr18 77398752 80373285 80089658 283627 2082861061
chr19 54075815 58617616 58440758 176858 1341657686
chr20 62148856 64444167 63944581 499586 1741460772
chr21 38164294 46709983 40088683 6621300 1203547746
chr22 36926516 50818468 39159777 11658691 988334122
awk 'BEGIN{OFS=FS="\t";}{if (NR>1){cov=$2/$4;print cov}}' count_all.txt >cov_percent.txt #计算覆盖率
awk 'BEGIN{OFS=FS="\t";}{if (NR>1){cov=$6/$2;print cov}}' count_all.txt >mean_depth.txt
sed -i '1 i cov_percent' cov_percent.txt
paste count_all.txt mean_depth.txt cov_percent.txt >count_all_done.txt
cat count_all_done.txt
Chrom coverage length Cut_N_length N_Number sum_of_depth mean_depth cov_percent
chrM 16569 16569 16568 1 34591025 2087.7 1.00006
chrX 151545600 156040895 154893130 1147765 4093132439 27.0092 0.978388
chrY 4004642 57227415 26415183 30812232 125469175 31.3309 0.151604
chr1 224300366 248956422 230481193 18475229 6482954037 28.903 0.973183
chr2 235466712 242193529 240548238 1645291 6376037101 27.0783 0.978875
chr3 194391189 198295559 198100139 195420 5286606057 27.1957 0.981277
chr4 185704231 190214555 189752667 461888 5330867456 28.7062 0.978665
chr5 176555626 181538259 181265378 272881 4949134249 28.0316 0.974017
chr6 164731880 170805979 170078724 727255 4420827669 26.8365 0.968563
chr7 155128428 159345973 158970132 375841 4173524521 26.9037 0.975834
chr8 141501464 145138636 144768136 370500 3842370943 27.1543 0.977435
chr9 118343356 138394717 121790590 16604127 3228529809 27.281 0.971695
chr10 129732215 133797422 133263135 534287 3626998764 27.9576 0.973504
chr11 131193733 135086622 134533742 552880 3548159595 27.0452 0.975173
chr12 130177621 133275309 133138039 137270 3490214776 26.8112 0.977764
chr13 95749291 114364328 97983128 16381200 2618593673 27.3484 0.977202
chr14 87847379 107043718 90568149 16475569 2359781279 26.8623 0.969959
chr15 81867181 101991189 84641348 17349841 2098133774 25.6285 0.967224
chr16 78578328 90338345 81805944 8532401 2341297030 29.7957 0.960545
chr17 78271424 83257441 82921074 336367 1911455487 24.4209 0.943927
chr18 77398752 80373285 80089658 283627 2082861061 26.9108 0.966401
chr19 54075815 58617616 58440758 176858 1341657686 24.8107 0.92531
chr20 62148856 64444167 63944581 499586 1741460772 28.0208 0.971917
chr21 38164294 46709983 40088683 6621300 1203547746 31.536 0.951997
chr22 36926516 50818468 39159777 11658691 988334122 26.7649 0.942971
知识补充：
1depth 测序深度
测序深度是指测序得到的总碱基数与待测基因组大小的比值，可以理解为基因组中每个碱基被测序到的平均次数。测序深度 = reads长度 × 比对的reads数目 / 参考序列长度。假设一个基因大小为2M，测序深度为10X，那么获得的总数据量为20M。
之所以测序深度=reads长度×比对的reads数目 / 参考序列长度，可以这么去理解，因为深度指的是每个碱基被测序到的平均次数。所以假设reads长度特别长，长到等于参考序列长度（当然这不可能）。那么，这时候reads数目就直接等于测序深度，但是reads长度远远没有参考序列长，所以才有了相当于reads长度/参考序列长度，这样的一个可以理解为矫正系数。
如果我们的研究目的仅仅是进行群体遗传分析等类型的分析。由于分析是通过计算基因组不同区域的多样性变化来推断进化选择压力，或构建进化树推断遗传分化关系，SNP检出率达到70~90%足以达到此类目的。所以，此类研究常见的测序深度选择在10X左右。
如果想检测个体的全基因组突变，来寻找某个特定功能突变，则我们推荐大于30X的测序深度，以保证接近于“毫无疏漏”。
如果测序样本是混样，例如，癌组织样本（癌细胞的细胞异质性很大），BSA分析中的群体混合池样本。那么在保证30X的测序深度的基础上，如果经费许可50X的测序深度更佳，以保证对混合池中低频变异的检测成功率（例如，新出现或即将消亡的癌细胞突变）。
2.coverage 覆盖度：
覆盖度是指测序获得的序列占整个基因组的比例。指的是基因组上至少被检测到1次的区域，占整个基因组的比例。当然，有些文章中也会将测序深度称为Coverage，容易给我们带来混淆。因此，我们还是需要根据语境来推断Coverage的意思。
我们假设基因组大小为G, 假定每次测序可从基因组任何位置上随机检测一个碱基。那么对于基因组上某一个固定碱基位置，在一次测序中，该位置被命中的概率为P （P=1/G）。我们将试验重复n次，相当于产生了n个碱基（n=c*G, c为测序深度）。 碱基的深度分布，相当于求该位置被测到0次，1次，…,n次的概率各是多少? 位点被检测到的n次的概率符合泊松分布。当然，由于概率极低，检测次数很大，所以这样的泊松分布接近于正态分布。
由于基因组中的高GC、重复序列等复杂结构的存在，测序最终拼接组装获得的序列往往无法覆盖有所的区域，这部分没有获得的区域就称为Gap。例如一个细菌基因组测序，覆盖度是98%，那么还有2%的序列区域是没有通过测序获得的。
上面统计所有染色体测序深度用的是python的脚本，这里也提供一个shell脚本，但是代码会相对长一点。
#!/bin/bash
for i in $(seq 1 22) #22对染色体
do
grep "chr${i} " /home/ubuntu/WGS_new/output/bwa_mapping/samtools_depth/KPGP-00001_L1.sorted.bam.depth |cut -f 3 |awk '{sum += $1;}END{print sum;}' >/home/ubuntu/WGS_new/output/bwa_mapping/samtools_depth/chr${i}.txt #注意END 的使用 |注意\t键的使用，脚本中可以直接键入tab，命令行下先Ctrl+V，再按tab键。这里之所以要加tab键，是因为，比如，如果>不加tab键，grep "chr1" 得到包括chr1，chr11，chr1_...>等。
sed -i "1 i chr${i}" /home/ubuntu/WGS_new/output/bwa_mapping/samtools_depth/chr${i}.txt # 1表示第一行，i表示插入，即在第一行前面插入新行，名称为chr${i}
done
grep "chrX " /home/ubuntu/WGS_new/output/bwa_mapping/samtools_depth/KPGP-00001_L1.sorted.bam.depth |cut -f 3 |awk '{sum += $1;}END{print sum;}' > /home/ubuntu/WGS_new/output/bwa_mapping/samtools_depth/chrx.txt #x染色体
grep "chrY " /home/ubuntu/WGS_new/output/bwa_mapping/samtools_depth/KPGP-00001_L1.sorted.bam.depth |cut -f 3 |awk '{sum += $1;}END{print sum;}' > /home/ubuntu/WGS_new/output/bwa_mapping/samtools_depth/chry.txt #y染色体
sed -i "1 i chrx" /home/ubuntu/WGS_new/output/bwa_mapping/samtools_depth/chrx.txt #添加列标题 chrx
sed -i "1 i chry" /home/ubuntu/WGS_new/output/bwa_mapping/samtools_depth/chry.txt
注意：awk命令要注意尽量冒号小写
awk '{sum += $1} END{print sum}' chr2.tmp
结果
3399
而
awk "{sum += $1} END{print sum}" chr2.tmp
报错
awk: cmd. line:1: {sum += } END{print sum}
awk: cmd. line:1: ^ syntax error
小总结
cut -f 可以提取特定列
cut -f 2 提取第2列
sed -n 2p 提取2行
sed -n 2,4p 提取2至3行
对于上诉得到的汇总结果存到count.txt中，并进行可视化，脚本如下
pdf()
library(ggplot2)
data <- read.table("count.txt",header = T)
#因为上一个命令，header=T已经将首行作为列标题，所以$1没法识别，用$加列名
data <-data[-1,] #将线粒体那行删除
data$Chrom <-factor(data$Chrom,levels = c(1:22,"X","Y")) #这一步很重要，这样横坐标就按顺序排列了
Chr <- data$Chrom
M <- data$mean_depth
C <- data$cov_percent
dt1 <- data.frame(Chr,M)
dt2 <- data.frame(Chr,C)
p1 <- ggplot(dt1,aes(x=Chr,y=M))+xlab("chrom ID")+ylab("mean_depth")+labs(title="KPGP-00001")+ylim(0,100)+geom_bar(aes(fill=Chr),stat = "identity")+theme_bw()+theme(plot.title = element_text(hjust = 0.5),panel.grid.major =element_blank(),panel.grid.minor = element_blank())
# 1 ggplot(dt1,aes(x=Chr,y=M,group=Chr))是构建一个画板
# 2 +xlab("chrom ID")+ylab("mean_depth") 是更改横纵坐标标签
# 3 +labs(title="KPGP-00001") 是增加主标题
# 4 ylim(0,100) 设置y轴刻度为1~100
# 5 +geom_bar(aes(fill=Chr),stat = "identity")是画出柱状图
# 6 theme_bw()为传统的白色背景和深灰色的网格线
# 7 theme是参数的设置 plot.title = element_text(hjust = 0.5) 是将title移到中间
# panel.grid.major =element_blank(),panel.grid.minor = element_blank()
#去除网格线
p2 <- ggplot(dt2,aes(x=Chr,y=C,group=1))+xlab("chrom ID")+ylab("coverage_percent")+labs(title="KPGP-00001")+ylim(0,1.5)+geom_line(aes(color=Chr))+geom_point(aes(color=Chr))+theme_bw()+theme(plot.title = element_text(hjust = 0.5),panel.grid.major =element_blank(),panel.grid.minor = element_blank())
p1
p2
dev.off()
最终结果如下
image.png
image.png
Y染色体的覆盖度特别低，不知道什么原因。
2.把bam文件按染色体分成小文件，并用IGV查看
bam文件太大了，有时候把bam文件拆分成小的文件是有必要的
这里用到一个工具bamtools,它的功能很多，这里选用split的功能
usage: bamtools [--help] COMMAND [ARGS]
Available bamtools commands:
convert Converts between BAM and a number of other formats
count Prints number of alignments in BAM file(s)
coverage Prints coverage statistics from the input BAM file
filter Filters BAM file(s) by user-specified criteria
header Prints BAM header information
index Generates index for BAM file
merge Merge multiple BAM files into single file
random Select random alignments from existing BAM file(s), intended more as a testing tool.
resolve Resolves paired-end reads (marking the IsProperPair flag as needed)
revert Removes duplicate marks and restores original base qualities
sort Sorts the BAM file according to some criteria
split Splits a BAM file on user-specified property, creating a new BAM output file for each value found
stats Prints some basic statistics from input BAM file(s)
See 'bamtools help COMMAND' for more information on a specific command.
这里只用lane1的bam文件做一个测试
nohup bamtools split -in KPGP-00001_L1.sorted.bam -reference &
#-reference指的是按reference拆分
拆分后，选取文件大小较小的chr19的bam文件，和chrY的文件（正好看看上一部布统计的覆盖度很低是不是统计错误）
#排序
samtools sort KPGP-00001_L1.sorted.REF_chrY.bam >chr_Y.sort.bam
samtools sort KPGP-00001_L1.sorted.REF_chr19.bam >chr_19.sort.bam
#建索引
samtools index chr_Y.sort.bam
samtools index chr_19.sort.bam
然后把文件导入IGV
先看一下chrY
image.png
可见覆盖度确实很低，并不是之前统计的问题
再看看chr19
image.png
chr19的覆盖程度明显大很多
3.bam格式转化为bw格式看测序深度分布
bam文件很大，很多时候才看很不方便。有些时候，我们只需要我们的测序数据在全基因组的测序深度这一个值，并不需要具体某条reads的碱基序列，碱基质量值。这样只需要把bam文件压缩为bw格式。
同样是只用lane1做测试
这里用到的一个软件deeptools
用到它其中的一个工具bamCoverage
#首先还是要先建立索引
nohup samtools index KPGP-00001_L1.sorted.bam &
再就是转换格式
nohup bamCoverage -b KPGP-00001_L1.sorted.bam -o KPGP-00001_L1.sorted.bw &
bw格式的大小才200M，可以轻松的导入IGV，可以从整体看出测序的深度和覆盖度情况。
比如
image.png
很显然，y染色体覆盖度很小
4.测序深度和GC含量的关系
首先，把全基因组的bam文件用 mpileup模式输出，根据 1000bp 的窗口滑动来统计每个窗口的测到的碱基数，GC碱基数，测序总深度。因为涉及到测序深度，所以先把各个lane的bam文件merge，再拿merge后的bam文件得到结果的前面一百万行的数据做测试。
nohup samtools merge KPGP-00001.sorted.merge.bam KPGP-00001_L*.sorted.bam &
samtools mpileup -f /home/ubuntu/WGS_new/input/genome/new/hg38.fa KPGP-00001.sorted.merge.bam KPGP-00001_L*.sorted.bam |head -n 1000000 >1M_mpileup.txt
先看一下输出文件的大概格式(这个文件是之前用lane2的bam做的测试结果)
image.png
包括6列：
分别是：
染色体号、碱基位置、参考序列对应的碱基、比对上的reads数目（也就是depth）、比对的情况、比对的质量
命令中之所以加一个-f 参数就是为了提供参考序列fasta，使最终得到的结果文件的第3列是对应的碱基，如果不加，则变成了N，如下
image.png
既然加了参考序列，就序列提前对参考序列构建索引，软件要求的是faidx index。
回到正题
现在统计测序深度和GC含量的关系，理论上测序得到的reads 拼接的结果应该是和参考序列一致的，因此统计参考序列的GC含量实际上就是统计对应的测序本身的insert 片段的GC含量。然后与此同时统计测序的对应区间的测序深度的值，观察两者的关系。
python脚本如下
import sys
args=sys.argv
filename=args[1]
count_depth={}
count_base={}
count_number={}
count_all={}
with open (filename) as fh:
for line in fh:
lineL=line.strip().split("\t")
chr_ID=lineL[0]
base_position = int(lineL[1])
window_position = int(base_position/1000)
base=lineL[2]
base.upper() #大写
depth=int(lineL[3])
chr_win=chr_ID+" "+str(window_position)
if chr_win not in count_base :
count_base[chr_win]=0
# if chr_win not in count_depth: 因为两个字典是按同一个key，所以这条可以省略
count_depth[chr_win]=0
count_number[chr_win]=0
if base=="G" or base=="C":
count_base[chr_win]+=1
count_depth[chr_win]+=depth
count_number[chr_win]+=1
count_all[chr_win]=str(count_base[chr_win])+' '+str(count_depth[chr_win])+' '+str(count_number[chr_win]) #把3个字典的value合并成一个字符串作为一个新的字典的value
for k,v in count_all.items():
print(k,v,sep="\t")
这个脚本写的很长，我相信肯定改进的地方很大，但是这是我自己写的脚本，所以凑活用吧。
#总结一下：python字典计数的方式可大致分为两种写法
# 1. if key not in dict:
# dict[key]=1 或者其他初始值
# else:
# dict[key]+=1
# 2. if key not in dict:
# dict[key]='' 或者 0
# dict[key]+=1
#两种用法虽然没什么大的差别，但是有的情况下，2能用的起来，而1不行或更麻烦，比如此脚本情况
# 其实区别在于1中的if else，执行了if就不会执行else。而2中执行了if会继续执行下面的语句
python3 count.py 1M_mpileup.txt |sort -k2,2n >count_all.txt
head count_all.txt
---------------
chr1 9 0 3 1
chr1 10 0 425756 529
chr1 11 133 35354 952
chr1 12 600 46231 1000
chr1 13 575 41108 1000
chr1 14 582 51064 1000
chr1 15 537 30959 1000
chr1 16 541 28821 1000
chr1 17 621 51877 1000
chr1 18 537 34340 1000
-------------------------------------
第一列是染色体号，
第二列是对应的窗口第几个窗口（从10开始）
第三列是对应GC碱基数
第四列是总测序深度
第五列是每个窗口统计到的碱基总数
再分别计算一下GC%和mean depth就可以可视化了
在这之前还是先看看刚统计的结果，可以发现，以1000bp的窗口滑动，统计到的碱基数目并不为1000，有的甚至很低，其实这很容易理解。我的理解是，对于参考基因组，仍然是存在一些含N的重复序列，而这列序列在比对的过程中是不计算的，samtools mpileup 生成的文本中这部分就已经去除了，所以碱基的位置信息不是连续的。而脚本中1000bp是按碱基的位置信息除以1000取int来做的，所以统计到的碱基数目不一致是正常的。
接下来给出可视化的R脚本：
library(ggplot2)
dt <- read.table("count.txt") #读入数据
#head(dt)
dt$GC <- dt$V3/dt$V5 #得到GC含量
dt$depth <-dt$V4/dt$V5 #得到平均的depth
dt <- dt[dt$GC>0,] #把GC等于0的部分去除
dt <- dt[dt$depth<100,] #把depth大于100的去除
lm_eqn <- function(x,y){
m <- lm(y ~ x);
eq <- substitute(italic(y) == a + b %.% italic(x)*","~~italic(r)^2~"="~r2,
list(a = format(coef(m)[1], digits = 2),
b = format(coef(m)[2], digits = 2),
r2 = format(summary(m)$r.squared, digits = 3)))
as.character(as.expression(eq));
} #添加线性回归方程
p <- ggplot(dt, aes(x = GC, y = depth)) +xlim(c(0,0.8))+ylim(c(0,100))+
geom_smooth(method = "lm", se=FALSE, color="black", formula = y ~ x)+
geom_point()+geom_text(x=0.6,y=80,label = lm_eqn(dt$GC , dt$depth), parse = TRUE)+theme_bw()+theme(panel.grid = element_blank())
#作图
结果图如下
image.png
从这个结果看似乎是负相关，所以再取其他百万行的结果重复一遍，看看是否得到同样的结果
这是第2000000至3000000之间的图
image.png
负相关，不知道什么原因，先不去管它
Part4 去重复（或标记重复序列）
以上对比对的横向探索结束
接下来继续纵向走，下一步是
去重复（或标记重复序列）
这里用两种软件：samtools和picard分别处理
1.samtools
这里本来应该是用samtools的markdup来标记或去除PCR重复，本来脚本已经写好开始运行了，但是报错了
脚本如下
#!/bin/bash
dir1=~/WGS_new/output/bwa_mapping
dir2=~/WGS_new/output/bwa_mapping/Dup
for i in $(seq 1 6)
do
samtools markdup -r $dir1/KPGP-00001_L${i}.sorted.bam $dir2/KPGP-00001_L${i}.sorted.rm.bam && echo "rm KPGP-00001_L${i} done"
samtools markdup -S $dir1/KPGP-00001_L${i}.sorted.bam $dir2/KPGP-00001_L${i}.sorted.mk.bam && echo "mk KPGP-00001_L${i} done"
done
报错信息如下
[markdup] error: no ms score tag.
[markdup] error: no ms score tag.
[markdup] error: no ms score tag.
[markdup] error: no ms score tag.
[markdup] error: no ms score tag.
[markdup] error: no ms score tag.
[markdup] error: no ms score tag.
仔细看了一下说明书
Usage: samtools markdup <input.bam> <output.bam>
Option:
-r Remove duplicate reads
-l INT Max read length (default 300 bases)
-S Mark supplemenary alignments of duplicates as duplicates (slower).
-s Report stats.
-T PREFIX Write temporary files to PREFIX.samtools.nnnn.nnnn.tmp.
-t Mark primary duplicates with the name of the original in a 'do' tag. Mainly for information and debugging.
--input-fmt-option OPT[=VAL]
Specify a single input file format option in the form
of OPTION or OPTION=VALUE
-O, --output-fmt FORMAT[,OPT[=VAL]]...
Specify output format (SAM, BAM, CRAM)
--output-fmt-option OPT[=VAL]
Specify a single output file format option in the form
of OPTION or OPTION=VALUE
--reference FILE
Reference sequence FASTA FILE [null]
-@, --threads INT
Number of additional threads to use [0]
The input file must be coordinate sorted and must have gone through fixmates with the mate scoring option on.
最后一行说到，必须gone through fixmates with the mate scoring option on
所以需要先fixmate
samtools fixmate KPGP-00001_L1.sorted.bam KPGP-00001_L1.sorted.fix.bam
然而又报错了
[bam_mating_core] ERROR: Coordinate sorted, require grouped/sorted by queryname
bam文件需要按quername排序，而samtools默认是按position排序的，而原始的bam文件为了减少储存，已经被我删除了。。。。
所以这一步没办法做了，因而samtools markdup去重就不做了，直接用picard的吧
2.picard
这里用GATK4的中的picard
#!/bin/bash
dir1=~/WGS_new/output/bwa_mapping
dir2=~/WGS_new/output/bwa_mapping/GATK/Dup
for i in $(seq 1 6)
do
#将重复序列删除，注意--REMOVE_DUPLICATES=true
gatk MarkDuplicates --REMOVE_DUPLICATES=true -I $dir1/KPGP-00001_L${i}.sorted.bam -M $dir2/KPGP-00001_L${i}.rm_metrics.txt -O $dir2/KPGP-00001_L${i}.rm.bam && echo "rm KPGP-00001_L${i} done"
#只是将重复标记，可以让后续软件识别，但是并没有删除
gatk MarkDuplicates -I $dir1/KPGP-00001_L${i}.sorted.bam -M $dir2/KPGP-00001_L${i}.mk_metrics.txt -O $dir2/KPGP-00001_L${i}.mk.bam && echo "mk KPGP-00001_L${i} done"
done
这里先解释为什么要去除重复，这部分节选自解螺旋的矿工公众号
首先，我们需要先理解什么是重复序列，它是如何产生的，以及为什么需要去除掉？要回答这几个问题，我们需要再次理解在建库和测序时到底发生了什么。
我们已经知道，在NGS测序之前都需要先构建测序文库：通过物理（超声）打断或者化学试剂（酶切）切断原始的DNA序列，然后选择特定长度范围的序列去进行PCR扩增并上机测序。
因此，这里重复序列的来源实际上就是由PCR过程中所引入的。因为所谓的PCR扩增就是把原来的一段DNA序列复制多次。可是为什么需要PCR扩增呢？如果没有扩增不就可以省略这一步了吗？
情况确实如此，但是很多时候我们构建测序文库时能用的细胞量并不会非常充足，而且在打断的步骤中也会引起部分DNA的降解，这两点会使整体或者局部的DNA浓度过低，这时如果直接从这个溶液中取样去测序就很可能漏掉原本基因组上的一些DNA片段，导致测序不全。而PCR扩增的作用就是为了把这些微弱的DNA多复制几倍乃至几十倍，以便增大它们在溶液中分布的密度，使得能够在取样时被获取到。所以这里大家需要记住一个重点，PCR扩增原本的目的是为了增大微弱DNA序列片段的密度，但由于整个反应都在一个试管中进行，因此其他一些密度并不低的DNA片段也会被同步放大，那么这时在取样去上机测序的时候，这些DNA片段就很可能会被重复取到相同的几条去进行测序。看到这里，你或许会觉得，那没必要去除不也应该可以吗？因为即便扩增了多几次，不也同样还是原来的那一段DNA吗？直接用来分析对结果也不会有影响啊！难道不是吗？
会有影响，而且有时影响会很大！最直接的后果就是同时增大了变异检测结果的假阴和假阳率。主要有几个原因：
DNA在打断的那一步会发生一些损失，主要表现是会引发一些碱基发生颠换变换（嘌呤-变嘧啶或者嘧啶变嘌呤），带来假的变异。PCR过程会扩大这个信号，导致最后的检测结果中混入了假的结果；
PCR反应过程中也会带来新的碱基错误。发生在前几轮的PCR扩增发生的错误会在后续的PCR过程中扩大，同样带来假的变异；
对于真实的变异，PCR反应可能会对包含某一个碱基的DNA模版扩增更加剧烈（这个现象称为PCR Bias）。因此，如果反应体系是对含有reference allele的模板扩增偏向强烈，那么变异碱基的信息会变小，从而会导致假阴。
PCR对真实的变异检测和个体的基因型判断都有不好的影响。GATK、Samtools、Platpus等这种利用贝叶斯原理的变异检测算法都是认为所用的序列数据都不是重复序列（即将它们和其他序列一视同仁地进行变异的判断，所以带来误导），因此必须要进行标记（去除）或者使用PCR-Free的测序方案（这种方案目前正变得越来越流行，特别是对于RNA-Seq来说尤为重要，现在著名的基因组学研究所——Broad Institute，基本都是使用PCR-Free的测序方案）。
跑完上面的流程的结果如下：
image.png
可以看出只是标记重复的（mk）文件是要比去除重复（rm）文件要大的。
再看一下之前的bam 文件
image.png
很显然这两个处理后的文件都是要比原始文件要大的。
分别看看mk和rm生成的文件有什么不同
samtools view KPGP-00001_L1.mk.bam |head
结果如下
image.png
samtools view KPGP-00001_L1.rm.bam |head
image.png
samtools view KPGP-00001_L1.sorted.bam |head
image.png
很明显，相比原始bam文件，rm多了一些markduplicate的字符，所以文件变大了
而，相比rm文件，mk，在flag值一列变大了，所以mk比rm文件大，这些通过文件内容很容易看出来
从下一步开始就直接用mk文件
因为如果 MarkDuplicates 这步把重复去掉，则会对 mate 信息产生影响，可以考虑在这步之后加上 FixMateInformation。
但是如果只是 MarkDuplicates而没有去除重复，则不会对mate信息产生影响，理论上不需要 FixMateInformation这一步
这个可以在哪得到论证呢
我们可以看看上面对3个bam文件查看的内容，对比flag值
（荧光标出来的部分）
mk.bam rm.bam sort.bam
1187 99 163
99 163 99
1123 163 99
163 99 163
1123 99 99
image.png
image.png
image.png
结果显而易见
Part5 碱基质量重校正
先下载一些必备数据库
#!/bin/bash
wget ftp://gsapubftp-anonymous@ftp.broadinstitute.org/bundle/hg38/1000G_phase1.snps.high_confidence.hg38.vcf.gz
wget ftp://gsapubftp-anonymous@ftp.broadinstitute.org/bundle/hg38/1000G_phase1.snps.high_confidence.hg38.vcf.gz.tbi
wget ftp://gsapubftp-anonymous@ftp.broadinstitute.org/bundle/hg38/Mills_and_1000G_gold_standard.indels.hg38.vcf.gz
wget ftp://gsapubftp-anonymous@ftp.broadinstitute.org/bundle/hg38/Mills_and_1000G_gold_standard.indels.hg38.vcf.gz.tbi
wget ftp://gsapubftp-anonymous@ftp.broadinstitute.org/bundle/hg38/hapmap_3.3.hg38.vcf.gz
wget ftp://gsapubftp-anonymous@ftp.broadinstitute.org/bundle/hg38/hapmap_3.3.hg38.vcf.gz.tbi
wget ftp://gsapubftp-anonymous@ftp.broadinstitute.org/bundle/hg38/dbsnp_146.hg38.vcf.gz
wget ftp://gsapubftp-anonymous@ftp.broadinstitute.org/bundle/hg38/dbsnp_146.hg38.vcf.gz.tbi
结果如下
total 5.1G
-rw-rw-r-- 1 ubuntu ubuntu 1.8G Aug 11 00:52 1000G_phase1.snps.high_confidence.hg38.vcf.gz
-rw-rw-r-- 1 ubuntu ubuntu 2.1M Aug 11 01:07 1000G_phase1.snps.high_confidence.hg38.vcf.gz.tbi
-rw-rw-r-- 1 ubuntu ubuntu 3.2G Aug 11 03:04 dbsnp_146.hg38.vcf.gz
-rw-rw-r-- 1 ubuntu ubuntu 2.4M Aug 11 03:20 dbsnp_146.hg38.vcf.gz.tbi
-rw-rw-r-- 1 ubuntu ubuntu 837 Aug 10 23:17 download.sh
-rw-rw-r-- 1 ubuntu ubuntu 60M Aug 11 01:08 hapmap_3.3.hg38.vcf.gz
-rw-rw-r-- 1 ubuntu ubuntu 1.5M Aug 11 01:08 hapmap_3.3.hg38.vcf.gz.tbi
-rw-rw-r-- 1 ubuntu ubuntu 20M Aug 11 01:07 Mills_and_1000G_gold_standard.indels.hg38.vcf.gz
-rw-rw-r-- 1 ubuntu ubuntu 1.5M Aug 11 01:07 Mills_and_1000G_gold_standard.indels.hg38.vcf.gz.tbi
在此之前先把之前下载的hg38.fa构建dict,注意dict和之前的bwa构建的index不一样
cd /home/ubuntu/WGS_new/input/genome/new
nohup gatk CreateSequenceDictionary -R hg38.fa -O hg38.dict &
创建了字典之后就可以继续了
#!/bin/bash
snp=/home/ubuntu/WGS_new/output/bwa_mapping/GATK/resource/dbsnp_146.hg38.vcf.gz
indel=/home/ubuntu/WGS_new/output/bwa_mapping/GATK/resource/Mills_and_1000G_gold_standard.indels.hg38.vcf.gz
#dir1=~/WGS_new/output/bwa_mapping
dir2=~/WGS_new/output/bwa_mapping/GATK/Dup
#dir3=/home/ubuntu/WGS_new/output/bwa_mapping/GATK/index
dir4=/home/ubuntu/WGS_new/output/bwa_mapping/GATK/resource
dir5=/home/ubuntu/WGS_new/output/bwa_mapping/GATK/recal
ref=~/WGS_new/input/genome/new/hg38.fa
for i in $(seq 1 6)
do
gatk BaseRecalibrator -R $ref -I $dir2/KPGP-00001_L${i}.mk.bam --known-sites $snp --known-sites $indel -O $dir5/KPGP-00001_L${i}.mk_recal.table
done
得到的table文件只是一个中间文件，还需要一步才能得到BQSR后的bam文件。
#!/bin/bash
snp=/home/ubuntu/WGS_new/output/bwa_mapping/GATK/resource/dbsnp_146.hg38.vcf.gz
indel=/home/ubuntu/WGS_new/output/bwa_mapping/GATK/resource/Mills_and_1000G_gold_standard.indels.hg38.vcf.gz
#dir1=~/WGS_new/output/bwa_mapping
dir2=~/WGS_new/output/bwa_mapping/GATK/Dup
#dir3=/home/ubuntu/WGS_new/output/bwa_mapping/GATK/index
dir4=/home/ubuntu/WGS_new/output/bwa_mapping/GATK/resource
dir5=/home/ubuntu/WGS_new/output/bwa_mapping/GATK/recal
ref=~/WGS_new/input/genome/new/hg38.fa
#for i in $(seq 1 6)
#do
#gatk BaseRecalibrator -R $ref -I $dir2/KPGP-00001_L${i}.mk.bam --known-sites $snp --known-sites $indel -O $dir5/KPGP-00001_L${i}.mk_recal.table
#done
for i in $(seq 1 6)
do
gatk ApplyBQSR -R $ref -I $dir2/KPGP-00001_L${i}.mk.bam -bqsr $dir5/KPGP-00001_L${i}.mk_recal.table -O $dir5/KPGP-00001_L${i}.mk_bqsr.bam
done
好了，到了这一步，利用生成的高质量的bam，接下来就是用4种工具来找变异，分别是GATK、bcftools、freebayes、varscan
再这之前先将各个lane生成的bam文件合并
nohup samtools merge KPGP-00001_bqsr.merge.bam *.bam
Part6 找变异
1.GATK
先对上一步生成的merge的bam文件建立索引
nohup samtools index KPGP-00001_bqsr.merge.bam &
然后是使用GATK的HaplotypeCaller
2.bcftools
#!/bin/bash
snp=/home/ubuntu/WGS_new/output/bwa_mapping/GATK/resource/dbsnp_146.hg38.vcf.gz
indel=/home/ubuntu/WGS_new/output/bwa_mapping/GATK/resource/Mills_and_1000G_gold_standard.indels.hg38.vcf.gz
#dir1=~/WGS_new/output/bwa_mapping
dir2=~/WGS_new/output/bwa_mapping/GATK/Dup
#dir3=/home/ubuntu/WGS_new/output/bwa_mapping/GATK/index
dir4=/home/ubuntu/WGS_new/output/bwa_mapping/GATK/resource
dir5=/home/ubuntu/WGS_new/output/bwa_mapping/GATK/recal
dir6=/home/ubuntu/WGS_new/output/bwa_mapping/GATK/vcf
ref=~/WGS_new/input/genome/new/hg38.fa
samtools mpileup -ugf $ref $dir5/KPGP-00001_bqsr.merge.bam | bcftools call --thread 6 -v -m -O z -o $dir6/KPGP-00001_bqsr.vcf.gz
