call variants from wheat RNA_seq

上周我们推送的何中虎研究员的文章（中国小麦产业发展与科技进步——小麦里我见过的最豪华作者阵容），目前阅读人数是8400多人，也是迄今为止我们公众号阅读人数最高的推送之一，这也直接让我们的关注人数猛增到3000人。过去一段时间，有不少小麦育种老师和专家也关注了我们，我们非常欢迎各位分享育种方面的经验。

今天要说一些说从RNAseq数据里得到序列变异的步骤。首先要交代一下背景，我们要研究一个EMS突变体（已回交多次），前期已经将突变基因定位到一个染色体区间，根据中国春参考序列，我们已经发现一个候选基因已经在水稻里被报道过了，测序发现，该基因确实发生了变异。后期安排了RNAseq实验，想在机制方面做一个有益的尝试。另外我们也想通过这样一个RNA_seq数据，比较在定位区间内还有那些基因发生了变异，因此就有了今天这样一个推送。

这里还要特别强调一点，这不是混池数据，分析过程中的一些参数请根据实验目的调整。

使用STAR将reads mapping至小麦基因组，然后使用sentieon流程（本质是GATK）call variant，接着使用SnpSift筛选高质量SNP，结合EMS诱变的特点，进一步排除可能的假阳性SNP，最后获得大概300个SNP，使用SnpEff注释SNP。根据前期遗传定位结果，我们发现只有4个SNP位于我们的区间内(20Mb)，但是只有一个SNP导致蛋白提前终止，该SNP所在的基因其在水稻里的直系同源基因已被报道，突变之后与我们的突变体表型非常相似。

下面是具体的流程，如果有兴趣欢迎交流。这个需要具有一定的高通量数据分析基础，其他的就没有什么特别的地方了。

# 工作目录
cd /data/rna_seq/genome/
#构建 A 基因组 index
STAR --runThreadN 10--runMode genomeGenerate --genomeDir ./--genomeFastaFiles CS_A_genome_part.fasta --sjdbFileChrStartEnd TGACv1_part_A.ss --limitGenomeGenerateRAM 68800833920

#!/usr/bin/env python
# -*- coding: utf-8 -*-
__author__ ='wheatomics'
import subprocess
# inpu.txt放着fastq文件的名字
with open('input.txt','r')as f:
    for line in f:
        line = line.strip().split()
        fq1, fq2 = line
        print fq1, fq2
        # 1. Mapping reads with STAR
        proc = subprocess.Popen(
            ['STAR','--twopassMode','Basic','--genomeDir','/data2/masw_data/rna_seq/dqyRNA-seq/masw_analysis/genome/','--runThreadN','20','--limitSjdbInsertNsj','5000000','--outSAMtype','BAM','SortedByCoordinate','--twopass1readsN','-1','--readFilesCommand','zcat','--sjdbOverhang','100','--outFilterMismatchNmax','10','--readFilesIn', fq1, fq2,
             '--outSAMattrRGline','ID:'+ fq1.split('_')[0],'SM:'+ fq1.split('_')[0],'PL:ILLUMINA'], shell=False)
        proc.wait()
        proc = subprocess.Popen(['mv','Aligned.sortedByCoord.out.bam','sorted1.bam'], shell=False)
        proc.wait()
        proc = subprocess.Popen(['samtools','view','-@','10','-h','-q','255','-o','sorted.bam','-O','BAM','sorted1.bam'], shell=False)
        proc.wait()
        proc = subprocess.Popen(['sentieon','util','index','sorted.bam'], shell=False)
        proc.wait()
        # 2. Metrics
        proc = subprocess.Popen(['sentieon','driver','-r','/data2/masw_data/rna_seq/dqyRNA-seq/masw_analysis/genome/CS_A_genome_part.fasta','-t','20','-i','sorted.bam','--algo','MeanQualityByCycle','mq_metrics.txt','--algo','QualDistribution','qd_metrics.txt','--algo','GCBias','--summary','gc_summary.txt','gc_metrics.txt','--algo','AlignmentStat','--adapter_seq',"''",'aln_metrics.txt','--algo','InsertSizeMetricAlgo','is_metrics.txt'], shell=False)
        proc.wait()
        proc = subprocess.Popen(['sentieon','plot','metrics','-o', fq1.split('_')[0]+'-metrics-report.pdf','gc=gc_metrics.txt','qd=qd_metrics.txt','mq=mq_metrics.txt','isize=is_metrics.txt'], shell=False)
        proc.wait()
        # 3. Remove Duplicate Reads
        proc = subprocess.Popen(['sentieon','driver','-t','20','-i','sorted.bam','--algo','LocusCollector','--fun','score_info','score.txt'], shell=False)
        proc.wait()
        proc = subprocess.Popen(['sentieon','driver','-t','20','-i','sorted.bam','--algo','Dedup','--rmdup','--score_info','score.txt','--metrics','dedup_metrics.txt','deduped.bam'], shell=False)
        proc.wait()
        # 4. Split reads at Junction
        proc = subprocess.Popen(['sentieon','driver','-r','/data2/masw_data/rna_seq/dqyRNA-seq/masw_analysis/genome/CS_A_genome_part.fasta','-t','20','-i','deduped.bam','--algo','RNASplitReadsAtJunction','--reassign_mapq','255:60','splitted.bam'], shell=False)
        proc.wait()
        # 5. Indel realigner
        proc = subprocess.Popen(['sentieon','driver','-r','/data2/masw_data/rna_seq/dqyRNA-seq/masw_analysis/genome/CS_A_genome_part.fasta','--read_filter','MapQualFilter,min_map_qual=40','-t','20','-i','splitted.bam','--algo','Realigner', fq1.split('_')[0]+'.realigned.bam'], shell=False)
        proc.wait()

Call SNP

#此处只统计了unique mapped的reads
sentieon driver -r /data2/masw_data/rna_seq/dqyRNA-seq/masw_analysis/genome/CS_A_genome_part.fasta --read_filter MapQualFilter,min_map_qual=60-t 10-i WT.realigned.bam -i br.realigned.bam --algo Genotyper--emit_conf 20--call_conf 20 WT_br_UG.vcf

筛选SNP

# 
# 要注意，EMS诱变的变异一般是C/T和A/G的变异，其他类型的变异频率很低很低。另外，突变里的关键功能变异理论上与参考序列不同。EMS mutations result in G:C to A:T mutations, whereas false positives could be any change. Thus, we retained only the alleles that corresponded to G:C to A:T mutations using SnpSift
cat WT_ms_UG.vcf | java -jar /data/snpEff/SnpSift.jar filter "(QUAL > 30) &(MQ > 40) & (QD > 5) & (FS < 30.0) & GEN[0].DP > 5 & GEN[1].DP > 5 & (((GEN[0].GT = '0/0') & (GEN[1].GT = '1/1')) | ((GEN[0].GT = '1/1') & (GEN[1].GT = '0/0')))"> WT_ms_UG_filtered.vcf
#上述筛选参数不是固定的，要根据实验和分析结果调整。具体每个参数表示的意思，请Google搜索SnpSift即可。

vcf文件里的一条染色体还是拆开的，需要合并成一个整体。

#!/usr/bin/env python
# -*- coding: utf-8 -*-
__author__ ='wheatomics'
chr =[['chr1A',471304005],['chr1B',438720154],['chr1D',452179604],['chr2A',462376173],['chr2B',453218924],
       ['chr2D',462216879],['chr3A',454103970],['chr3B',448155269],['chr3D',476235359],['chr4A',452555092],
       ['chr4B',451014251],['chr4D',451004620],['chr5A',453230519],['chr5B',451372872],['chr5D',451901030],
       ['chr6A',452440856],['chr6B',452077197],['chr6D',450509124],['chr7A',450046986],['chr7B',453822637],
       ['chr7D',453812268]]
with open('WT_br_UG_filtered.vcf','r')as f:
    for line in f:
        if line.startswith('#'):
            print line,
        else:
            line = line.replace('_part1','')
            line = line.strip().split('')
            if line[0].endswith('part2'):
                for i in chr:
                    if line[0].split('_')[0]== i[0]:
                        line[1]= int(line[1])+ int(i[1])
                        line[0]= line[0].split('_')[0]
            for m in line[:-1]:
                print str(m)+'',
            print line[-1]+'',

最好再修改下vcf的表头信息

##contig=<ID=chr1A,length=594102056,assembly=unknown>
##contig=<ID=chr2A,length=780798557,assembly=unknown>
##contig=<ID=chr3A,length=750843639,assembly=unknown>
##contig=<ID=chr4A,length=744588157,assembly=unknown>
##contig=<ID=chr5A,length=709773743,assembly=unknown>
##contig=<ID=chr6A,length=618079260,assembly=unknown>
##contig=<ID=chr7A,length=736706236,assembly=unknown>
##contig=<ID=chrUn,length=480980714,assembly=unknown>

转换之后就要统计SNP的信息，比如染色体上的SNP个数等

#使用snpEff注释SNP
java -Xmx8g-jar snpEff.jar IWGSCv1.0 WT_itr_UG_filtered_whole.vcf > WT_itr_UG_filtered_whole_eff.vcf

根据上述结果，即可进行下一步的分析。

SnpEff这个我们前面也介绍过，可以参考使用SnpEff 对SNP结果进行分析。

欢迎关注“小麦研究联盟”，了解小麦新进展