这两篇文章所用材料都有小麦微核心种质，如果结合两个结果是不是有意思的结果出来呢？实际上第一篇文章也是可以做关联分析的。下面我们就一步一步来，实际上这也是我的实验记录。

首先要在NCBI上下载数据，这一步略去不表，参见我们以前的推送SRA数据的下载以及在线blast--或许与你了解的不一样这第二步就是将下载的数据mapping至中国春基因组，并获得包含变异信息的gvcf文件。因为是转录组数据，所以mapping软件使用的是STAR，这里没有列出如何将基因组序列index。下面是流程，整个从流程从sra文件开始，到bam文件结束。具体应用到自己的项目上时，要根据需要修改。还是那句话，小麦里的变异分析，我也是新媳妇上花轿——头一回。

#!/usr/bin/env python
# -*- coding: utf-8 -*-
__author__ ='wheatomics'
import subprocess
with open('input.txt','r')as f:
    for line in f:
        line = line.strip().split()
        sra, rg = line
        print sra, rg
        proc = subprocess.Popen(['fastq-dump','--split-3','--defline-qual','+','--defline-seq','@$ac-$si/$ri','--helicos',  sra +'.sra'], shell=False)
        proc.wait()
        proc = subprocess.Popen(
                ['fastp','-w','10','-i', sra +'_1.fastq','-o', sra +'_out.1.fastq','-I', sra +'_2.fastq',
                 '-O', sra +'_out.2.fastq'], shell=False)
        proc.wait()
        # 1. Mapping reads with STAR
        proc = subprocess.Popen(
            ['STAR','--twopassMode','Basic','--genomeDir','/data2/Fshare/FastaAndIndex/IWGSC_v1.0_STAR/',
             '--runThreadN','20','--limitSjdbInsertNsj','5000000','--outSAMtype','BAM','SortedByCoordinate','--twopass1readsN','-1',
             '--sjdbOverhang','100','--outFilterMismatchNmax','6','--readFilesIn', sra +'_out.1.fastq', sra +'_out.2.fastq',
             '--outSAMattrRGline','ID:'+ rg,'SM:'+ rg,'PL:ILLUMINA'], shell=False)
        proc.wait()
        proc = subprocess.Popen(['mv','Aligned.sortedByCoord.out.bam','sorted.bam'], shell=False)
        proc.wait()
        proc = subprocess.Popen(['shred','-u','-z', sra +'_1.fastq'], shell=False)
        proc.wait()
        proc = subprocess.Popen(['shred','-u','-z', sra +'_2.fastq'], shell=False)
        proc.wait()
        proc = subprocess.Popen(['shred','-u','-z', sra +'_out.1.fastq'], shell=False)
        proc.wait()
        proc = subprocess.Popen(['shred','-u','-z', sra +'_out.2.fastq'], 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/Fshare/FastaAndIndex/IWGSC_v1.0_STAR/IWGSC_v1.0_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', rg +'-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/Fshare/FastaAndIndex/IWGSC_v1.0_STAR/IWGSC_v1.0_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/Fshare/FastaAndIndex/IWGSC_v1.0_STAR/IWGSC_v1.0_part.fasta',
                                 '-t','20','-i','splitted.bam','--algo','Realigner', rg +'realigned.bam'], shell=False)
        proc.wait()

得到bam文件之后，通过下面的命令就可以获得vcf文件

# 此处为命令示例，多个bam文件，可使用 -i 添加
sentieon driver -r /data2/Fshare/FastaAndIndex/IWGSC_v1.0_STAR/IWGSC_v1.0_part.fasta --read_filter MapQualFilter,min_map_qual=60-t 10-i Aimengniu.realigned.bam -i Aodesa3.realigned.bam -i Baibiansui.realigned.bam -i Baidatou.realigned.bam 

下面这一步就是筛选SNP了，筛选这一步要到下期介绍了。后面会多一些发表文章常见的结果。因为本文是90个材料的RNA_seq的数据，所以会包括群体遗传分析和进化方面的一些内容，下面推送会包括以为内容，GWAS，eGWAS，结合其他已知信息尝试锁定QTL候选基因等。

最后祝大家假期愉快，万事如意！还有，情人节快乐，单身的早日找到另一半