Ms1基因组重测序数据分析

Hi, 大家好！我是胖丫，今天的推送让我来。这篇很早就写好了，Ms1文章发布后，我们专门组织了一个再分析群，借此来学习和构建一个基本的分析流程，为以后开展相关课题做个基础。群里还有另外一个小伙伴广伟也在尝试这样的一个工作，他所使用的流程与我的不一样，原本想着一起推送出来，但终究没有如愿，人家早早的跑出结婚了，借此也恭喜他们。小麦基因组巨大，下载测序数据就有大概600G，分析过程中产生的一系列中间文件也是巨大无比，整个过程持续了一月之久，好在是计算机干活。本次分析重点在于打通分析流程，精力有限，所有程序的参数或者筛选条件等等均未优化过，所以在使用过程中大家还需根据结果多探索。本文也只是列出了代码以及简单的注释，没有过多的去讨论和说明很多额外的东西，主要还是因为我也是一个二把刀，也是初次尝试分析小麦重测序数据，也没有啥经验。所以，文中若有不足之处，万望指出。

言归正传，我们前面推送了三篇与Ms1有关的文章，不了解的小伙伴可以点击下面的链接了解。

PNAS在线发表MS1的图位克隆工作--方法值得借鉴

兼听则明，有容乃大--- 由最近Ms1的文章想到的

call variations from wheat RNA_seq --以MS1为例

大家也看到了，前面我们已经分析过RNA_seq数据了。由于我们已经知道了是哪个基因的突变以及突变位点，所以我们首先查看了ms1e突变体和野生型的Ms1（TraesCS4B01G017900.1）基因的序列，然而，我们并未发现报道的突变位点有变异，再看了一些突变位点的reads数目，发现较少。所以，我们推测这可能是因为覆盖度不够导致的。所以就得看全基因组重测序数据了。

首先需要下载sra文件，将不含后缀的文件名（即不含.sra）保存在一个名为“dnainput.txt”的文件中,每行一个sra文件名。本脚本需要的程序有：

1、pfastq-dump

2、BWA （BWA index网上教程多的是，未列出）

3、sentieon

4、运行上述软件所需要的背景软件

输入文件： dnainput.txt

SRR5873698
SRR5873699
SRR5873700
SRR5873701
SRR5873702
SRR5873703
SRR5873704
SRR5873705

输出文件： .g.vcf.gz .g.vcf.gz.tbi

提醒：本文只是建立一个基本的流程，具体到每个项目中，还需要根据目的调整和优化相关流程或参数。本代码一次只mapping一个样本的数据，这里列出的是野生型的mapping过程，最终生成g.vcf文件。

#!/usr/bin/env python
# -*- coding: utf-8 -*-
__author__ = 'wheatomics'
importsubprocess
withopen('dna_input.txt', 'r') asf:
    sra_input = ''
    ref_fasta = '/data/IWGSC_v1.0_bwa/161010_Chinese_Spring_v1.0_pseudomolecules_parts.fasta'
    sample = 'wild'  # sample name
    forlineinf:
        sra = line.strip().split()[0]
        printsra+'.sra'
        # sra to fastq.gz
        proc = subprocess.Popen(['pfastq-dump', '--threads', '10', '--outdir', './', sra+'.sra', '--split-3', '--gzip'], shell=False)
        proc.wait()
        # Mapping each set of input fastq with BWA-MEM, sorting
        proc = subprocess.Popen('bwa mem -M -R \'@RG\tID:'+sra+'\tSM:'+sample+'\tPL:ILLUMINA\''+
                                 ' -t 12 -K 10000000 '+ref_fasta+' '+sra.split('.')[0] +'_1.fastq.gz'+' '+sra.split('.')[0] +'_2.fastq.gz'+
                                 ' | '+'sentieon util sort -r '+ref_fasta+' -o '+sra+'.sorted.bam'+' -t 10 --sam2bam -i -', shell=True)
        proc.wait()
        # delete all the fastq.gz
        proc = subprocess.Popen(['shred', '-u', '-z', sra.split('.')[0] +'_1.fastq.gz'], shell=False)
        proc.wait()
        proc = subprocess.Popen(['shred', '-u', '-z', sra.split('.')[0] +'_2.fastq.gz'], shell=False)
        proc.wait()
        # get all bam file and put them into next step
        sra_input = sra_input+' -i '+sra+'.sorted.bam'
        proc.wait()
    printsra_input
    # Remove Duplicate Reads on the multiple sorted BAM files
    proc = subprocess.Popen('sentieon driver -t 10'+sra_input+' --algo LocusCollector --fun score_info score.txt', shell=True)
    proc.wait()
    proc = subprocess.Popen('sentieon driver -t 10'+sra_input+  ' --algo Dedup --rmdup --score_info score.txt --metrics dedup_metrics.txt deduped.bam', shell=True)
    proc.wait()
    # Indel realigner
    proc = subprocess.Popen(['sentieon', 'driver', '-r', ref_fasta, '-t', '10', '-i', 'deduped.bam', '--algo', 'Realigner',
                             'realigned.bam'], shell=False)
    proc.wait()
    # Base recalibration
    proc = subprocess.Popen(['shred', '-u', '-z', 'deduped.bam'], shell=False)
    proc.wait()
    proc = subprocess.Popen(['sentieon', 'driver', '-r', ref_fasta, '-t', '10', '-i', 'realigned.bam', '--algo',
                             'QualCal', 'recal_data.table'], shell=False)
    proc.wait()
    proc = subprocess.Popen(['sentieon', 'driver', '-r', ref_fasta, '-t', '10', '-i', 'realigned.bam', '-q',
                             'recal_data.table', '--algo', 'QualCal', 'recal_data.table.post'], shell=False)
    proc.wait()
    proc = subprocess.Popen(['sentieon', 'driver', '-t', '10', '--algo', 'QualCal', '--plot', '--before',
                             'recal_data.table', '--after', 'recal_data.table.post', 'recal.csv'], shell= False)
    proc.wait()
    proc = subprocess.Popen(['sentieon', 'plot', 'bqsr', '-o', 'recal_plots.pdf', 'recal.csv'], shell=False)
    proc.wait()
    # HC Variant caller
    proc = subprocess.Popen(['sentieon', 'driver', '-r', ref_fasta, '-t', '10', '-i', 'realigned.bam', '-q',
                             'recal_data.table', '--algo', 'Haplotyper', '--emit_mode', 'gvcf', '--emit_conf',
                             '10', '--call_conf', '30', sample+'-hc.g.vcf.gz'], shell=False)
    proc.wait()
    proc = subprocess.Popen(['shred', '-u', '-z', 'realigned.bam'], shell=False)
    proc.wait()

合并多个gvcf生成vcf文件;

sentieon driver -r /data/IWGSC_v1.0_STAR/IWGSC_v1.0_part.fasta -t10--algo GVCFtyper Ms1_DNA.vcf mutant-hc.g.vcf.gz wild-hc.g.vcf.gz

因为我们使用的参考序列染色体被分成了两部分，这里要先合并起来。

python combine_vcf.py > Ms1_DNA_whole.vcf

### combine_vcf.py
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]]
withopen('Ms1_DNA.vcf', 'r') asf:
    forlineinf:
        ifline.startswith('#'):
            printline,
        else:
            line = line.replace('_part1', '')
            line = line.strip().split('\t')
            ifline[0].endswith('part2'):
                foriinchr:
                    ifline[0].split('_')[0] == i[0]:
                        line[1] = int(line[1]) +int(i[1])
                        line[0] = line[0].split('_')[0]
            forminline[:-1]:
                printstr(m) +'\t',
            printline[-1] +'\n',

下面一步就是筛选SNP，用来计算突变体和野生型的SNP index以及两者的比值。筛选这一步我也没有啥经验，参考文章标准的标准以及自己摸索了一些标准，感觉都不是太理想，代码也就不列出了，精力有限，也就没在继续折腾下去。好在，从结果来看，还是可以看到4B上有峰存在，也即能够检测到突变位点位于4B染色体的某一段，如下图。

至此，分析告一段落，下面说一下有意思的发现。首先，RNA_seq能检测到的多态，但是在DNA_seq里却没有，猜测在RNA反转录过程中会引入一些错误。第二就是RNA_seq测序寻找候选基因基因时一定要选对时间点。第三是关于Ms1这个基因的，如下图，此处能检测到ms1e的突变，除了这些，同时还发现，在Ms1基因下游的表达量要远高于上游，貌似是紧挨着Ms1的一个基因，该基因并未注释出来，通常认为与Ms1是在同一个转录本上。另外，还发现Ms1下游也有降解组reads支持，这说明此处也许有small RNA参与调控，上述表达丰度的不同可能就是此处的small RNA导致的。最终是不是与雄性不育有关还需要更多的实验证据验证。期待未来某一天有人能够阐明此处的生物学作用。