GATK做变异检测

一、测序数据的准备

新建工作目录，GATK-analysis

mkdir RNASeq-analysis
cd RNASeq-analysis/

新建data文件夹；
mkdir data

1.1 公司返回的测序下机数据

例如：6个NGS测序样本；
Raw data: Sample1_R1.fq.gz; Sample1_R2.fq.gz; Sample2_R1.fq.gz; Sample2_R2.fq.gz; Sample3_R1.fq.gz; Sample3_R2.fq.gz; Sample4_R1.fq.gz; Sample4_R2.fq.gz; Sample5_R1.fq.gz; Sample5_R2.fq.gz; Sample6_R1.fq.gz; Sample6_R2.fq.gz;

1.2 从SRA数据库下载数据

第一步，获取SRA编号Accession: PRJNA/SRP；

第二步，下载安装sra-tookit

conda install sra-tools

第三步，下载数据

prefetch SRR******

第四步，对数据进行拆分

使用 SRA-toolkit 软件包里的 fastq-dump 对数据进行拆分，并转换为 fastq 格式。
nohup fastq-dump --split-3 SRR****** 1>SRR******.log 2>SRR******.err &
将测序数据全部置于data文件夹，
mv Sample* data/

二、对数据进行质控和过滤

for file in `ls *_R1.fq.gz | perl -lpe 's/_R1.fq.gz//'`; do fastp -i ${file}_R1.fq.gz -I ${file}_R2.fq.gz -o ${file}_R1.clean.fq.gz -O ${file}_R2.clean.fq.gz; done

三、参考基因组准备

主要准备基因组文件Species.genome.fasta;
新建references文件夹：
mkdir references
将参考文件置于references文件夹；
mv Species.genome.fasta ~/GATK-analysis/references/

四、bwa比对

新建bwa工作目录：

 mkdir bwa
cd bwa '
vim bwa.sh

#bwa.sh
bwa index ../references/Species.genome.fasta
for file in `ls ../data/*_R1.fq.gz | perl -lpe 's/_R1.fq.gz//'`; do 
    bwa mem -t 4 -R '@RG\tID:foo\tPL:illumina\tSM:${file}' ../references/Species.genome.fasta ${file}_R1.fq.gz ${file}_R2.fq.gz > ${file}.sam
    samtools view -b ${file}.sam > ${file}.bam
    samtools sort ${file}.bam > ${file}.sorted.bam
    samtools index ${file}.sorted.bam
done

运行bwa:
nohup bash bwa.sh &

五、GATK calling SNP

直接下载安装：
wget https://github.com/broadinstitute/gatk/releases/download/4.1.7.0/gatk-4.1.7.0.zip
欲了解详情，参见官方文档：
GATK 4.1.7
新建工作目录：
mkdir gatk
cd gatk
将bwa生成的排序文件移至gatk目录，并建立参考基因组的软连接：

mv ../bwa/*.sorted.bam .
ln -s ../references/Species.genome.fasta .

编写gatk脚本：
vim gatk.sh

samtools faidx Species.genome.fasta
gatk CreateSequenceDictionary -R Species.genome.fasta -O Species.genome.dict
for file in *.sorted.bam;do
    gatk MarkDuplicates -I ${file} -O ${file}.markdup.bam -M ${file}.markdup_metrics.txt
    samtools index ${file}.markdup.bam
    gatk HaplotypeCaller -R Species.genome.fasta --emit-ref-confidence GVCF -I ${file}.markdup.bam -O ${file}.g.vcf
    gatk GenotypeGVCFs -R Species.genome.fasta -V ${file}.g.vcf -O ${file}.vcf    
    bgzip -f ${file}.vcf
    tabix -p vcf ${file}.vcf.gz
    gatk SelectVariants -select-type SNP -V ${file}.vcf.gz -O ${file}.snp.vcf.gz    
    gatk VariantFiltration -V ${file}.snp.vcf.gz --filter-expression "QD < 2.0 || MQ < 40.0 || FS > 60.0 || SOR > 3.0 || MQRankSum < -12.5 || ReadPosRankSum < -8.0" --filter-name "PASS" -O ${file}.snp.filter.vcf.gz
    gatk SelectVariants -select-type INDEL -V ${file}.vcf.gz -O ${file}.indel.vcf.gz
    gatk VariantFiltration -V ${file}.indel.vcf.gz --filter-expression "QD < 2.0 || FS > 200.0 || SOR > 10.0 || MQRankSum < -12.5 || ReadPosRankSum < -8.0" --filter-name "PASS" -O ${file}.indel.filter.vcf.gz
    gatk MergeVcfs -I ${file}.snp.filter.vcf.gz -I ${file}.indel.filter.vcf.gz -O ${file}.filter.vcf.gz
done

运行GATK：
nohup bash gatk.sh &

得到结果文件${file}.filter.vcf.gz，里面包含SNP和indel。

参考：

Shifu Chen, Yanqing Zhou, Yaru Chen, Jia Gu; fastp: an ultra-fast all-in-one FASTQ preprocessor, Bioinformatics, Volume 34, Issue 17, 1 September 2018, Pages i884–i890, https://doi.org/10.1093/bioinformatics/bty560
Li H., Handsaker B., Wysoker A., Fennell T., Ruan J., Homer N., Marth G., Abecasis G., Durbin R. and 1000 Genome Project Data Processing Subgroup (2009) The Sequence alignment/map (SAM) format and SAMtools. Bioinformatics, 25, 2078-9. [PMID: 19505943]
Li H A statistical framework for SNP calling, mutation discovery, association mapping and population genetical parameter estimation from sequencing data. Bioinformatics. 2011 Nov 1;27(21):2987-93. Epub 2011 Sep 8. [PMID: 21903627]
Genome Analysis Toolkit

最后编辑于：2020.05.16 23:12:18