在Windows下完成QTL-seq&MutMap

最新推荐文章于 2024-05-21 08:36:35 发布

Nuvolar

最新推荐文章于 2024-05-21 08:36:35 发布

阅读量1.5k

点赞数 2

文章标签： windows python 生物信息学

本文链接：https://blog.csdn.net/sinat_41244077/article/details/120811545

版权

纯粹就是记录下如何在Window下少写代码完成QTL-seq&MutMap,并没有详细讲解QTL-seq&MutMap的原理和每一步背后的含义。详细原理以及实验设计可以看文章
QTL-seq: rapid mapping of quantitative trait loci in rice by whole genome resequencing of DNA from two bulked populations.;
Genome sequencing reveals agronomically important loci in rice using MutMap
High performance pipeline for MutMap and QTL-seq.
QTL-seq 分析原理及流程
 使用MutMap快速定位突变基因：原理篇
QTL-seq和MutMap都是从BSA发展过来的，已经比较成熟，测序成本也不大，只需要两三个混合的样本(bulk)就可以了。简单地说就是把数量性状转换为质量性状，利用detla SNP这个指标来定位QTL。
本文主要是利用Iwate Biotechnology Research Center的pipelines。但在windows下运行会有一点点问题，需要修改下。

工具准备

conda（需要添加到环境变量, 建议Python版本在3.8之后，miniconda 就可以了）Windows 10下安装Miniconda3
perl5 (需要添加到环境变量, 我装的草莓版的)windows环境安装Perl
java (需要添加到环境变量)JAVA(windows)安装教程
SRA Toolkit (用于将sra文件转为reads)下载地址
bowtie2 (需要添加到环境变量，最新版的还没被编译到windows平台，这里用的是2.3.4版本，名字含mingw的zip就可以。运行bowtie2需要perl在环境变量中)下载地址
snpEff (用来注释vcf文件，需要java) 下载地址
Trimmomatic (用来去reads的接头，需要java)下载地址
samtools (需要添加到环境变量，这个是1.3.1版本。原作者Li Heng 写过能用mingw/msy2编译的makefile文件，但在最新版的1.13版的编译中总是报错，所以找了个旧的支持windows的版本)下载地址
QTL-seq 下载地址 (建议先下载到本地再用conda安装。在code 下download zip。) 这个库本身的功能就是call snp 和计算delta snp。但是这里用了subprocess和mutiprocess，所以有些地方需要按照windows环境进行修改，就是该QTL-seq文件夹下qtlseq的py文件。
1、trim.py和qtlplot.py由于Trimmomatic 和snpEff 需要用java -jar 命令来启动所以相关命令要进行修改。
trim.py 修改44行和65行命令 java -jar 后面跟trimmomatic-0.39.jar的具体存放位置就可以了。这个命令说明如果要去接头的话，只能用双端测序的结果
trim.py
qtlplot.py 修改36行。java -jar 后面跟snpEff.jar的具体存放位置。

2、由于本身的pipelines用的是bwa进行比对。试了一下发现，bwa在windwos下运行太慢了，一个30x的样本比对要27个小时，这里给它换成bowtie2来比对，同样的样本只要5个小时就可以完成了。这里需要修改refindex.py和alignment.py。
refindex 就是改了下构建引索的命令
替换的refindex.py

import sys
import subprocess as sbp
from qtlseq.utils import time_stamp, clean_cmd, call_log

class RefIndex(object):

    def __init__(self, args):
        self.args = args
        self.out = args.out

    def run(self):
        print(time_stamp(),
              'start to index reference fasta.',
              flush=True)

        cmd1 = 'bowtie2-build {
   0} {
   1}/10_ref/ref\
                >> {
   1}/log/bowtie2.log \
                2>&1'.format(self.args.ref, self.out)

        cmd2 = 'samtools faidx {
   0} \
                >> {
   1}/log/samtools.log \
                2>&1'.format(self.args.ref, self.out)

        cmd1 = clean_cmd(cmd1)
        cmd2 = clean_cmd(cmd2)

        try:
            sbp.run(cmd1,
                    stdout=sbp.DEVNULL,
                    stderr=sbp.DEVNULL,
                    shell=True,
                    check=True)
        except sbp.CalledProcessError:
            call_log(self.out, 'bowtie2', cmd1)
            sys.exit(1)

        try:
            sbp.run(cmd2,
                    stdout=sbp.DEVNULL,
                    stderr=sbp.DEVNULL,
                    shell=True,
                    check=True)
        except sbp.CalledProcessError:
            call_log(self.out, 'samtools', cmd1)
            sys.exit(1)

        print(time_stamp(),
              'indexing of reference successfully finished.',
              flush=True)

alignment改了下比对命令，毕竟家用的windows不太可能有100g的运行内存。
替换的alignment.py

import sys
import os
import subprocess as sbp
from qtlseq.utils import time_stamp, clean_cmd, call_log


class Alignment(object):

    def __init__(self, args):
        self.args = args
        self.out = args.out

    def run(self, fastq1, fastq2, index):
        print(time_stamp(),
              'start to align reads of {} by bowtie2.'.format(index),
              flush=True)

        cmd1 = 'bowtie2 -p {0} -x {3}/10_ref/ref -1 {1} -2 {2} -S {3}/20_bam/{4}.sam >> {3}/log/alignment.log '.format(self.args.threads,fastq1,fastq2,self.out,index)
        cmd2 = 'samtools fixmate -O bam {0}/20_bam/{1}.sam {0}/20_bam/{1}_fixmate.bam '.format(self.out,index)
        cmd3 = 'rm {0}/20_bam/{1}.sam '.format(self.out,index)
        cmd4 = 'samtools sort -O bam -m 1G -@ {0} -o {1}/20_bam/{2}_sort.bam {1}/20_bam/{2}_fixmate.bam '.format(self.args.threads,self.out,index)
        cmd5 = 'rm {0}/20_bam/{1}_fixmate.bam '.format(self.out,index)
        cmd6 = 'samtools rmdup -sS {0}/20_bam/{1}_sort.bam {0}/20_bam/{1}_rmdup.bam '.format(self.out,index)
        cmd7 = 'rm {0}/20_bam/{1}_sort.bam '.format(self.out,index)
        cmd8 = 'samtools view -b -f 2 -F 2048 -o {0}/20_bam/{1}.bam {0}/20_bam/{1}_rmdup.bam '.format(self.out,index)
        cmd9 = 'rm {0}/20_bam/{1}_rmdup.bam '.format(self.out,index)

        cmd = [cmd1, cmd2, cmd3, cmd4, cmd5, cmd6, cmd7, cmd8, cmd9]
        for i in cmd:
            tem = clean_cmd(i)
            if tem.startswith("rm"):
                target = tem.split(" ")[1]
                os.remove(target)
                continue
            try:
                sbp.run(tem,
                    stdout=sbp.DEVNULL,
                    stderr=sbp.DEVNULL,
                    shell=True,
                    check=True)
            except sbp.CalledProcessError:
                call_log(self.out, 'alignment', tem)
                sys.exit(1)

        print(time_stamp(),
              'alignment {} successfully finished.'.format(index),
              flush=True)

3、在mpileup.py和vcf2index.py中用了pool.map 可能会有些问题，打不开子进程。还有管道符"|"可能在subprocess中用不了。这里改了下。
替换的mpileup.py

import sys
import re
import os
import glob
import subprocess as sbp
import multiprocessing
import multiprocessing.pool
from multiprocessing import Pool
from multiprocessing.dummy import Pool as ThreadPool
from qtlseq.vcf2index import Vcf2Index
from qtlseq.utils import time_stamp, clean_cmd, call_log

class Mpileup(object):
    def __init__(self, args):
        self.out = args.out
        self.args = args

    def get_bams(self, label):
        bams = glob.glob('{}/20_bam/{}*.bam'.format(self.out, label))
        return bams

    def merge(self):
        for label in ['parent', 'bulk1', 'bulk2']:
            bams = self.get_bams(label)
            if len(bams) == 1:
                path_to_bam = os.path.abspath(bams[0])
                cmd1 = 'ln -s {0} {1}/20_bam/{2}.unsorted.bam'.format(path_to_bam,
                                                                   self.out,
                                                                   label)
            else:
                cmd1 = 'samtools merge -f {
   0}/20_bam/{
   1}.unsorted.bam \
                                          {
   0}/20_bam/{
   1}*.bam \
                                          >> {
   0}/log/samtools.log \
                                          2>&1'.format(self.out, label)

            cmd2 = 'samtools sort -m {
   0} \
                                  -@ {
   1} \
                                  -o {
   2}/20_bam/{
   3}.bam \
                                  {
   2}/20_bam/{
   3}.unsorted.bam \
                                  >> {
   2}/log/samtools.log \
                                  2>&1'.format(self.args.mem,
                                               self.args.threads,
                                               self.out,
                                               label)

            cmd3 = 'samtools index {
   0}/20_bam/{
   1}.bam \
                                   >> {
   0}/log/samtools.log \
                                   2>&1'.format(self.out, label)

            cmd4 = 'rm {}/20_bam/{}.*.bam'.format(self.out, label)

            cmd = [cmd1, cmd2, cmd3, cmd4]
            for i in cmd:
                tem = clean_cmd(i)
                if tem.startswith("rm"):
                    target = glob.glob(tem.split(" ")[1])
                    for t in target:
                        os.remove(t)
                    continue
                if tem.startswith("ln"):
                    target = tem.split(" ")
                    os.symlink(target[2], target[3])
                    continue
                try:
                    sbp.run(tem,
                        stdout=sbp.DEVNULL,
                        stderr=sbp.DEVNULL,
                        shell=True,
                        check=True)
                except sbp.CalledProcessError:
                    call_log(self.out, 'samtools', tem)
                    sys.exit(1)

    def get_header(self):
        ref = open(self.args.ref, 'r')
        pattern = re.compile('>')
        chr_names = []
        for line in ref:
            if pattern.match(line):
                line = line.rstrip('\n')
                chr_name = re.split('[> ]',line)[1]
                chr_names.append(chr_name)
        return chr_names

    def mpileup(self, chr_name):
        print("call variants in {} ".format(chr_name))
        cmd1 = 'samtools mpileup -t AD,ADF,ADR -B -q {0} -Q {1} -C {2} -v -u -o {3}/30_vcf/tem_mpileup.{4}.vcf -r {4} -f {5} --ignore-RG {3}/20_bam/parent.bam {3}/20_bam/bulk1.bam {3}/20_bam/bulk2.bam >> {3}/log/bcftools.{4}.log 2>&1 '.format(self.args.min_MQ,self.args.min_BQ,self.args.adjust_MQ,self.out,chr_name,self.args.ref)
        cmd2 = 'bcftools call -vm -f GQ,GP -O z -o {1}/30_vcf/tem_call.{2}.vcf.gz {1}/30_vcf/tem_mpileup.{2}.vcf >> {1}/log/bcftools.{2}.log 2>&1 '.format(self.args.threads,self.out, chr_name) 
        cmd3 = 'bcftools filter -i "INFO/MQ>={3}" -O z -o {1}/30_vcf/qtlseq.{2}.vcf.gz {1}/30_vcf/tem_call.{2}.vcf.gz >> {1}/log/bcftools.{2}.log 2>&1 '.format(self.args.threads,self.out, chr_name, self.args.min_MQ) 
        cmd4 = 'tabix -f -p vcf {0}/30_vcf/qtlseq.{1}.vcf.gz >> {0}/log/tabix.{1}.log 2>&1 '.format(self.out, chr_name)
        cmd5 = 'rm {0}/30_vcf/tem_mpileup.{1}.vcf'.format(self.out, chr_name)
        cmd6 = 'rm {0}/30_vcf/tem_call.{1}.vcf.gz'.format(self.out, chr_name)

        cmd = [cmd1, cmd2, cmd3, cmd4, cmd5, cmd6]

        for i in cmd:
            tem = clean_cmd(i)
            if tem.startswith("rm"):
                target = glob.glob(tem.split(" ")[1])
                for t in target:
                    os.remove(t)
                continue
            try:
                sbp.run(tem,
                    stdout=sbp.DEVNULL,
                    stderr=sbp.DEVNULL,
                    shell=True,
                    check=True)
            except sbp.CalledProcessError:
                if tem.startswith("tabix"):
                    call_log(self.out, 'tabix.{0}'.format(chr_name), tem)
                if tem.startswith("bcftools") or tem.startswith("samtools"):
                    call_log(self.out, 'bcftools.{0}'.format(chr_name), tem)
                sys.exit(1)

    def concat(self):
        cmd1 = 'cat {0}/log/bcftools.*.log > {0}/log/bcftools.log'.format(self.out)
        cmd2 = 'cat {0}/log/tabix.*.log > {0}/log/tabix.log'.format(self.out)

        cmd3 = 'bcftools concat --threads {
   1} -O z \
                                -o {
   0}/30_vcf/qtlseq.vcf.gz \
                                {
   0}/30_vcf/qtlseq.*.vcf.gz \
                                >> {
   0}/log/bcftools.log \