一、测序数据的准备

新建工作目录，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。

---

参考：

1. 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
2. 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]
3. 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]
4. Genome Analysis Toolkit