使用工具GATK4。 GATK基础 RNA-seq分析完整版BSR参考 GATK官方RNA-seq calling流程最新版的BSR分析流程流程参考 WT 3个单株，混池。三个技术重复。 NS (实验组) 3个单株，混池。三个技术重复。 3个数量有点少，就暂且练习BSR吧。小麦的BSR BSR和BSA的比对方式不一致。基于DNA水平的重测序，可以测到所有的碱基变化情况，需要整个基因组比对。基于RNA表达水平的BSR,更多的是转录本的比对。相对而言，BSR的测序价格相对较低。

如果不懂参数的意义，不要随便修改参数，对于进程，尽量按照我的设置，有些程序你给的内存或进程太多，硬件没有这么高的配置有可能会报错。特别是基于Java的程序，很容易内存溢出。

下机数据： 1.质控

##获取输入的数组的长度
bef=(sample1 sample2 ... sampleN)
sample_num=${#bef[@]}

for ((i=0;i<$sample_num;i++));
do
inA1=../origin/${bef[$i]}$bg1;
inA2=../origin/${bef[$i]}$bg2;
out1=${bef[$i]}"paired-R1.fq.gz";
out2=${bef[$i]}"paired-R2.fq.gz";
unpaired1=${bef[$i]}"unpaired-R1.fq.gz";
unpaired2=${bef[$i]}"unpaired-R2.fq.gz";
fastp --thread 16 -i $inA1 -o $out1 -I $inA2 -O $out2 --html ${bef[$i]}".html" --json ${bef[$i]}".json" 2>${bef[$i]}fastp_out
echo success end  QC_${bef[$i]} fastp at `date`;
done

genome=/home/chaim/disk/maize/Zea_mays.B73_RefGen_v4.42.fa
gene_gtf=/home/chaim/disk/maize/Zea_mays.B73_RefGen_v4.42.gtf

比对到参考基因组STAR 2-pass-align 比对

2.比对（目前推荐的RNA比对工具是HISAT2和STAR） HISAT2的使用，在我的RNA-seq栏中有详细信息。此处使用STAR进行比对（STAR也被推荐用于RNA-seq call snp的比对）。 STAR的最新版下载安装下载安装，配置环境变量。 STAR github

STAR的安装

尽量使用最新版本的软件，因为旧版本可能会有莫名其妙的问题，同时新版支持使用STAR-fusion检测基因融合

#当前最新版本是2.7.3a，到上述链接github找最新版本
wget https://github.com/alexdobin/STAR/archive/2.7.3a.tar.gz
tar -xzf 2.7.3a.tar.gz
cd STAR-2.7.3a
cd bin/Linux_x86_64/
./STAR

配置环境变量 vim ~/.bash_profile 添加如下export PATH=$HOME/disk/STAR/STAR-2.7.3a/bin/Linux_x86_64:$PATH 重庆终端，STAR --version即可查看当前版本信息。 2.1构建基因组索引文件。（无论是hisat2还是STAR此步骤超慢）参考2-pass-mode 参考

STAR 参数介绍：

--runThreadN 6 使用线程数

--outSAMtype BAM SortedByCoordinate 将结果SAM文件排序并生成BAM，使用的软件和samtools效果类似(目前版本一定不要使用该参数，否则会因为cpu占用率过高，而程序挂掉，后续使用samtools完成这步)

--genomeDir ./starhuman150 上面建立的index目录

--readFilesIn AA_R1_trimmed.fastq AA_R2_trimmed.fastq 要输入的测序数据

--outFileNamePrefix AA 可以自定义生成文件的前缀

2.2 使用STAR进行2步比对

genome=/home/chaim/disk/maize/Zea_mays.B73_RefGen_v4.42.fa
gene_gtf=/home/chaim/disk/maize/Zea_mays.B73_RefGen_v4.42.gtf
bef=(sample1 sample2 ……samplen)
sample_num=${#bef[@]}
#构建基因组的索引
STAR --runThreadN 30 \
        --runMode genomeGenerate \
        --genomeDir ./stat-1-index \
        --genomeFastaFiles $genome \
        --sjdbGTFfile $gene_gtf \
        --sjdbOverhang 149 \
        2 >log

#在star>2.4.1版本后，最新版命令已经支持直接的两步比对。
for((i=0;i<sample_num;i++));
do
star_index=./stat-1-index
fq1=../origin/cleandata/${bef[$i]}"clean-R1.fq.gz"
fq2=../origin/cleandata/${bef[$i]}"clean-R2.fq.gz"
STAR --runThreadN 24 --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  ${bef[$i]}_
#自定义输出文件的前缀
echo star  end satr ${bef[$i]}`date`

2.3 sam转bam，排序构建索引基因组文件>500M,最好不要使用picard等java程序构建索引，否则后续GATK可能会报错。

for ((i=0;i<sample_num;i++));
do
mv ${bef[$i]}_Aligned.out.sam ${bef[$i]}.sam
#因为star输出排序后的bam文件时，系统cpu占用率太高，会报错。所以此处手动sort sam排序为bam
samtools sort -@ 24 -O bam -o ${bef[$i]}.bam ${bef[$i]}.sam
#生成bai文件，用于导入IGV
samtools index -@ 24 ${bef[$i]}".bam" ${bef[$i]}.bai
#统计比对的结果
samtools flagstat -@ 24 ${bef[$i]}.bam >${bef[$i]}.flagstat
done

生成bam文件后的处理去除重复的PCR可以使用的工具有picard,samtools,sambamba ,三个工具对比关于去除重复PCR的去除与否一般是RNA-seq不需要去除，call snp和chip-seq,dap-seq需要去除。原理详解使用picard去除重复的reads. java -Xms48g -jar $picard MarkDuplicates REMOVE_DUPLICATES= true INPUT=${bef[[$i]}.map.bam OUTPUT=${bef[[$i]}.repeatmark.bam METRICS_FILE=${bef[[$i]}.bam.metrics

3.0 合并同一组数据的bam文件（同一个样本的2次下机数据，最开始如果没有合并，此时应该合并）

java -Xmx48g -jar $picard MergeSamFiles I=2447-20_L1.sort.bam I=2447-20_L2.sort.bam O=2447-20.bam

<<!

截取部分数据用于测试数据

samtools view -h Z58.bam |head -20000 >tmp.sam
samtools sort -O bam -@ 4 tmp.sam -o tmp.bam
samtools index tmp.bam

#gatk IndexFeatureFile --feature-file /public/home/tang/chaim/maize/zea_mays.vcf
for((i=0;i<sample_num;i++));
do
samtools sort -O bam -@ 8 ${bef[$i]}.bam -o ${bef[$i]}.sort.bam
samtools index ${bef[$i]}.sort.bam
#添加头文件
java -jar $picard AddOrReplaceReadGroups I=${bef[$i]}.sort.bam O=${bef[$i]}.add.bam ID=${bef[$i]} LB=${bef[$i]} PL=illumina PU=${bef[$i]} SM=${bef[$i]} 
java -jar $picard MarkDuplicates I=${bef[$i]}.add.bam O=${bef[$i]}.rmd.bam CREATE_INDEX=true VALIDATION_STRINGENCY=SILENT M=${bef[$i]}.metrics 
#RNA-seq call snp 使用GATK增加的一步
gatk SplitNCigarReads -R $genome -I ${bef[$i]}.rmd.bam -O ${bef[$i]}.split.bam
gatk BaseRecalibrator -R $genome -I ${bef[$i]}.split.bam --use-original-qualities -O ${bef[$i]}.recal.table --known-sites $dbsnp_vcf
gatk ApplyBQSR -R $genome -I ${bef[$i]}.split.bam -bqsr ${bef[$i]}.recal.table -O ${bef[$i]}.recal.bam
gatk --java-options "-Xmx8g -Xms8g" HaplotypeCaller -R $genome -I ${bef[$i]}.recal.bam --dbsnp $dbsnp_vcf --native-pair-hmm-threads 8 -stand-call-conf 20 -O ${bef[$i]}.g.vcf.gz -ERC GVCF --verbosity ERROR
done
!

3.1 先处理有生物学重复的样品

bef=(sample1 sample2 sample3 ... sampleN)
sample_num=${#bef[@]}
for((i=0;i<$sample_num;i++));
do
#先合并转录组的生物学重复合bam文件
java -Xmx48g -jar $picard MergeSamFiles I=${bef[$i]}"-1.bam" I=${bef[$i]}"-2.bam" I=${bef[$i]}"-3.bam" O=${bef[$i]}.bam
#创建一个新的环境bioconda
#conda create -n bioconda
#用conda安装sambamba
#conda install -y -c bioconda sambamba 
#去除PCR重复
sample=${bef[$i]}
sambamba markdup --overflow-list-size 600000  --tmpdir='./' -r $sample.bam ${sample}_rmd.bam
#前面的STAR比对的时候没有添加头文件,此处增加Read group
gatk AddOrReplaceReadGroups -I ${sample}_rmd.bam -O ${sample}_rmd_add.bam -LB $sample -PL ILLUMINA -PU $sample -SM $sample
#创建索引，然后SplitNCigarReads
#gatk CreateSequenceDictionary -R $genome
gatk SplitNCigarReads -R $genome -I ${sample}_rmd_add.bam -O ${sample}_rmd_split.bam

#BQSR （Base Quality Recalibration）
gatk --java-options "-Xmx12G" BaseRecalibrator -I ${sample}_rmd_split.bam -R $genome -O ${sample}_recal.table --known-sites $dbsnp_vcf
gatk --java-options "-Xmx12G" ApplyBQSR -I ${sample}_rmd_split.bam -R $genome -O ${sample}_recal.bam -bqsr ${sample}_recal.table

#使用gvcf模式Variant Calling
bam=${sample}_recal.bam
gatk --java-options "-Xmx8G" HaplotypeCaller -ERC GVCF -R $genome -I $bam --dbsnp $dbsnp_vcf -O ${sample}.g.vcf
done

4.1再处理没有生物学重复的样品

pbef=(group1 group2 ... groupn)
sample_n=${#pbef[@]}
for((i=0;i<$sample_n;i++));
do
#创建一个新的环境bioconda
#conda create -n bioconda
#用conda安装sambamba
#conda install -y -c bioconda sambamba 
#去除PCR重复
sample=${pbef[$i]}
sambamba markdup --overflow-list-size 600000  --tmpdir='./' -r $sample.bam ${sample}_rmd.bam
#前面的STAR比对的时候没有添加头文件,此处增加Read group
gatk AddOrReplaceReadGroups -I ${sample}_rmd.bam -O ${sample}_rmd_add.bam -LB $sample -PL ILLUMINA -PU $sample -SM $sample
#创建索引，然后SplitNCigarReads
#gatk CreateSequenceDictionary -R $genome
gatk SplitNCigarReads -R $genome -I ${sample}_rmd_add.bam -O ${sample}_rmd_split.bam

#BQSR （Base Quality Recalibration）
gatk --java-options "-Xmx12G" BaseRecalibrator -I ${sample}_rmd_split.bam -R $genome -O ${sample}_recal.table --known-sites $dbsnp_vcf
gatk --java-options "-Xmx12G" ApplyBQSR -I ${sample}_rmd_split.bam -R $genome -O ${sample}_recal.bam -bqsr ${sample}_recal.table

#使用gvcf模式Variant Calling
bam=${sample}_recal.bam
gatk --java-options "-Xmx8G" HaplotypeCaller -ERC GVCF -R $genome -I $bam --dbsnp $dbsnp_vcf -O ${sample}.g.vcf
done

5.合并前面的gvcf到一个文件

gatk --java-options "-Xmx8G" CombineGVCFs -R $genome --dbsnp $dbsnp_vcf --variant sample1.g.vcf.gz --variant sample2.g.vcf.gz --variant group1.g.vcf.gz --variant group2.g.vcf.gz -O sample_group.g.vcf.gz

6 genotypeGVCFs

gatk --java-options "-Xmx12G" GenotypeGVCFs -R $genome -V sample_group.g.vcf.gz -O sample_group.common.vcf.gz 
##当有多个样本池的时候，在vcf_gop中添加多个sample_group名称即可。
vcf_gop=(sample_group1 ... sample_groupN)
group_n=${#pbef[@]}
for((i=0;i<$group_n;i++))
do
#7.filter variant
gatk --java-options "Xmx4G" VariantFiltration -V ${vcf_gop[$i]}.common.vcf.gz --window 35 --cluster 3 --filter-name "lowQD" --filter-expression "QD < 2.0" --filter-name "highFS" --filter-expression "FS > 30.0" -O ${vcf_gop[$i]}.Filter.vcf
#7.1select SNP
gatk --java-options "-Xmx12G" SelectVariants -R $genome -O ${vcf_gop[$i]}.filter_SNP.vcf --variant ${vcf_gop[$i]}.Filter.vcf --select-type-to-include SNP  2>${vcf_gop[$i]}.select_snp.err
#7.2 select indel
gatk --java-options "-Xmx12G" SelectVariants -R $genome -O ${vcf_gop[$i]}.filter_INDEL.vcf --variant ${vcf_gop[$i]}.Filter.vcf --select-type-to-include INDEL  2>${vcf_gop[$i]}.select_indel.err

#7.3 merge INDEL and SNP
#gatk MergeVcfs -I ${vcf_gop[$i]}.Filter_INDEL.vcf -I ${vcf_gop[$i]}.Filter_SNP.vcf -O ${vcf_gop[$i]}.Filter.vcf
#8. 转换vcf为table格式
gatk VariantsToTable -R $genome -V ${vcf_gop[$i]}.Filter.vcf -F CHROM -F POS -F REF -F ALT -GF AD -GF DP -GF GQ -GF PL -O ${vcf_gop[$i]}.table
done

9.变异注释

注释工具有vep,annova,GATK,snpEff

方法一：

#SnpEff version SnpEff 4.3T (build 2017-11-24 10:18), by Pablo Cingolani

snpEff="/public/home/tang/chaim/soft/snpEff/snpEff/snpEff.jar"
java -jar $snpEff ann zm442 sample_group.Filter.vcf -t -stats sample_group.summary.html -csvStats sample_group.csv >sample_group.ann.vcf

方法二：

#VariantQC是DISCVRSeq其中的一个函数，其是基于gatk的一个Java包。
#java -jar /public/home/tang/chaim/soft/DISCVRSeq/DISCVRSeq.jar VariantQC -O output_file -R ref.fa -V variant.vcf
DISCVRSeq=/public/home/tang/chaim/soft/DISCVRSeq/DISCVRSeq.jar
genome=/home/chaim/disk/maize/Zea_mays.B73_RefGen_v4.42.fa
#创建基因组文件的索引信息
samtools faidx $genome
java -jar $picard CreateSequenceDictionary R=$genome O=zm442.dict

#创建变异文件的tbi索引
#gatk输出vcf文件时会自动创建索引文件，一般情况不需要自己创建索引
#gatk IndexFeatureFile -F variant.vcf

#开始分析变异
##此处的maxContigs根据不同物种，尽量设置大的值，大于该物种的参考基因组中未组装到染色体上的contig数量。否则会报错。
java -jar $DISCVRSeq VariantQC -O group_sample.html -R $genome --maxContigs 52 -V group_sample.Filter.vcf

拿到变异文件后，后续和BSA分析一样。

最新工具速递----可变剪切贝叶斯分析不同的剪接BANDITS 基因融合检测工具STAR-fusion 参考文章 STAR-fusion的使用 WGCNA的数据挖掘多样本的WGCNA WGCNA官网教程

变异注释的方式参考地址

annovar下载地址 annovar需要教育邮箱，解决办法有很多。 snpEFF安装使用 snpEFF注释 java -Xmx128g -jar /home/snpEff.jar zm442 common_filtration.vcf > common.eff.vcf 会输出三个文件， snpEff_genes.txt snpEff_summary.html common.eff.vcf gatk的Funcotator目前只可用于人的变异注释 vcp用于注释变异安装需要root权限，依赖比较多。 VariantQC用于变异注释 VariantQC的产出是html文件，类似于FastQC。可视化图表比较美观。

写在最后

BSA和BSR分析结果不一致的问题。原因详细讲解我的BSA和BSR的结果并不一致。可能原因是： 1. 剪切式比对导致的SNP检测错误 2. RNA编辑的影响 3. 基因组印记现象的影响（表观遗传） 4. 目标基因不表达或低表达，以及无法检测非SNP的变异