一、測序數(shù)據(jù)的準(zhǔn)備
新建工作目錄针贬，GATK-analysis
mkdir RNASeq-analysis
cd RNASeq-analysis/
新建data文件夾；
mkdir data
1.1 公司返回的測序下機(jī)數(shù)據(jù)
例如：6個(gè)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數(shù)據(jù)庫下載數(shù)據(jù)
第一步逊朽，獲取SRA編號Accession: PRJNA******/SRP******；
第二步曲伊，下載安裝sra-tookit
conda install sra-tools
第三步叽讳，下載數(shù)據(jù)
prefetch SRR******
第四步追他，對數(shù)據(jù)進(jìn)行拆分
使用 SRA-toolkit 軟件包里的 fastq-dump 對數(shù)據(jù)進(jìn)行拆分，并轉(zhuǎn)換為 fastq 格式岛蚤。
nohup fastq-dump --split-3 SRR****** 1>SRR******.log 2>SRR******.err &
將測序數(shù)據(jù)全部置于data文件夾邑狸，
mv Sample* data/
二、對數(shù)據(jù)進(jìn)行質(zhì)控和過濾
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
三涤妒、參考基因組準(zhǔn)備
主要準(zhǔn)備基因組文件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
運(yùn)行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
運(yùn)行GATK：
nohup bash gatk.sh &
得到結(jié)果文件${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
