GATK4流程学习之RNA-Seq variant callin
2021-02-14 本文已影响0人
小贝学生信
GATK4流程学习之背景知识与前期准备 - 简书
GATK4流程学习之DNA-Seq variant calling(Germline：SNP+INDEL) - 简书
GATK4流程学习之RNA-Seq variant calling(SNP+INDEL) - 简书
补：Mutect2+scRNAseq+cancer cell - 简书
参考https://gatk.broadinstitute.org/hc/en-us/articles/360035531192-RNAseq-short-variant-discovery-SNPs-Indels-流程，以SRR11618713、SRR11618726、SRR11618727三个单细胞转录组数据为例，进行RNA-Seq variant calling(SNP+INDEL)分析；
由于基本流程类似DNA-seq，所以每一步详细介绍可参考上一篇笔记。
GATK:Best Practices Workflows
conda activate GATK
cd GATK/rna-seq
1、ascp下载fastq.gz数据
https://www.ebi.ac.uk/ena/browser/view/SRR11618713
https://www.ebi.ac.uk/ena/browser/view/SRR11618726
https://www.ebi.ac.uk/ena/browser/view/SRR11618727
mkdir -p ./raw/qc
cd raw
cat > ascp.link
fasp.sra.ebi.ac.uk:/vol1/fastq/SRR116/013/SRR11618713/SRR11618713_1.fastq.gz
fasp.sra.ebi.ac.uk:/vol1/fastq/SRR116/013/SRR11618713/SRR11618713_2.fastq.gz
fasp.sra.ebi.ac.uk:/vol1/fastq/SRR116/026/SRR11618726/SRR11618726_1.fastq.gz
fasp.sra.ebi.ac.uk:/vol1/fastq/SRR116/026/SRR11618726/SRR11618726_2.fastq.gz
fasp.sra.ebi.ac.uk:/vol1/fastq/SRR116/027/SRR11618727/SRR11618727_1.fastq.gz
fasp.sra.ebi.ac.uk:/vol1/fastq/SRR116/027/SRR11618727/SRR11618727_2.fastq.gz
cat ascp.link |while read sample
do
ascp -QT -l 300m -P33001 \
-i ~/miniconda3/envs/GATK/etc/asperaweb_id_dsa.openssh \
era-fasp@$sample .
done
2、质控过滤(optional)
fastqc *.gz -o ./qc
- 由于fastqc报告结果总体合格，未进行过滤。
- 利用trimmomatic过滤低质量fastq片段代码参考如下
trimmomatic PE -phred33 -trimlog logfile sample_1.fastq.gz sample_2.fastq.gz out.read_1.fq.gz \
out.trim.read_1.fq.gz out.read_2.fq.gz out.trim.read_2.fq.gz \
ILLUMINACLIP:/home/shensuo/biosoft/Trimmomatic/Trimmomatic-0.36/adapters/TruSeq3-PE.fa:2:30:10 \
SLsampleINGWINDOW:5:20 LEADING:5 TRAILING:5 MINLEN:50
3、star比对
mkdir ../bam/star_log ; cd ../bam
star_index=/home/data/gma49/GATK/ref/hg38/star-index/GRCh38_gencode_v22_CTAT_lib_Apr032020.plug-n-play/ctat_genome_lib_build_dir/ref_genome.fa.star.idx
cat > SRR.list
SRR11618713
SRR11618726
SRR11618727
cat SRR.list | while read sample
do
cd star_log
fq1=../../raw/${sample}_1.fastq.gz
fq2=../../raw/${sample}_2.fastq.gz
STAR --runThreadN 4 --genomeDir $star_index \
--twopassMode Basic --outReadsUnmapped None --chimSegmentMin 12 \
--alignIntronMax 100000 --chimSegmentReadGapMax parameter 3 \
--alignSJstitchMismatchNmax 5 -1 5 5 \
--readFilesCommand zcat --readFilesIn $fq1 $fq2 --outFileNamePrefix ${sample}_
cd ../ ;mv ./star_log/${sample}_Aligned.out.sam ./${sample}.sam
samtools sort -o ${sample}.bam ${sample}.sam
samtools index $sample.bam
samtools flagstat $sample.bam > $sample.flagstat
rm ${sample}.sam
done
1
4、去重mark duplicate
- 意义同上一篇DNA-seq笔记。这里使用的是其它教程推荐的sambamba软件，而不是之前用的gatk合并picard里的
MarkDuplicates功能. - 一开始使用sambamba的最新0.6.8版本，部分bam文件会出现
Segmentation fault (core dumped)报错，重新安装了0.6.6版本均可正常操作。
cat SRR.list | while read sample
do
sambamba markdup --overflow-list-size 10000 --tmpdir='./' -r $sample.bam ${sample}_rmd.bam
done
5、SplitNCigarReads
-
这一步是RNA-seq特异性的一步。因为mRNA转录本是主要由DNA的外显子exon可变剪切组合而成
2
- 把转录信息比对到基因组上时，比对到intron的区域均为N，直接根据此结果会出现假阳性variant结果。
- 所以需要将跨越了n个intron的read片段，切割为n+1个子read片段。
cat SRR.list | while read sample
do
gatk SplitNCigarReads -R $GENOME \
-I ./${sample}_rmd.bam \
-O ${sample}_rmd_split.bam
gatk AddOrReplaceReadGroups -I ${sample}_rmd_split.bam \
-O ${sample}_rmd_split_add.bam -LB $sample -PL ILLUMINA -PU $sample -SM $sample
done
AddOrReplaceReadGroups与之前bwa比对时的-R参数，添加必要的group信息
6、校正碱基质量分数 Base Quality Recalibration
cat SRR.list | while read sample
do
gatk BaseRecalibrator \
-I ${sample}_rmd_split_add.bam \
-R $GENOME \
--known-sites /home/data/gma49/GATK/ref/hg38/var_ref/dbsnp_146.hg38.vcf.gz \
--known-sites /home/data/gma49/GATK/ref/hg38/var_ref/Mills_and_1000G_gold_standard.indels.hg38.vcf.gz \
--known-sites /home/data/gma49/GATK/ref/hg38/var_ref/1000G_phase1.snps.high_confidence.hg38.vcf.gz \
-O ${sample}_recal.table
gatk ApplyBQSR \
--bqsr-recal-file ${sample}_recal.table \
-R /home/data/gma49/GATK/ref/hg38/star-index/GRCh38_gencode_v22_CTAT_lib_Apr032020.plug-n-play/ctat_genome_lib_build_dir/ref_genome.fa \
-I ${sample}_rmd_split_add.bam \
-O ${sample}_recal.bam
done
6、Variant Calling
- 为每个样本生成 g.vcf
mkdir ../vcf
cat SRR.list | while read sample
do
gatk HaplotypeCaller \
--native-pair-hmm-threads 10 \
-R $GENOME \
-D /home/data/gma49/GATK/ref/hg38/var_ref/dbsnp_146.hg38.vcf.gz \
-I ${sample}_recal.bam \
--minimum-mapping-quality 30 \
-ERC GVCF \
-O ../vcf/${sample}.g.vcf
done
- 合并 g.vcf
cd ../vcf
gatk CombineGVCFs \
-R $GENOME \
--variant SRR11618713.g.vcf \
--variant SRR11618726.g.vcf \
--variant SRR11618727.g.vcf \
-O SRR3.all.g.vcf
gatk GenotypeGVCFs \
-R $GENOME \
-V SRR3.all.g.vcf \
-O SRR3.vcf
至此完成了初步的RNA-Seq variant calling(SNP+INDEL)流程分析，接下来可对vcf结果进行进一步过滤，其中也会涉及到很多知识点。会另外整理总结。
