细菌全基因组测序数据构建SNP系统发育树

2019-06-19 本文已影响0人涤生生

本次分析使用到的数据及数据处理流程参考文献“https://doi.org/10.1186/s13756-018-0352-y ”。考虑到时间及存储有限，且此次分析的目的仅为跑pipeline，因此本次测试仅用了ERP020566 中的小部分数据，见ERR.list。

整个数据处理过程包括如下步骤：

数据下载及转化：包括NCBI下载sra数据以及候选的参考基因组数据；数据预处理
参考基因组选择
call SNP
构建系统发育树

吐槽文章的一个错误。文章中call snp后是把基因组中的噬菌体序列及重复序列区域去掉的（After excluding repetitive regions and mobile genetic elements, 159,154 SNP positions (length of the reference genome: 5,306,618 bp) ），而在图注中标明159,154 SNP 构建的该菌核心基因组系统发育树（The image shows a ML tree based on 159.154 SNPs of the core genome. A K. variicola isolate ）。这个说法是有问题的，在查阅资料后，使用了一套新的方法构建核心基因组snp系统发育树，新方法会另写文章。

1. 数据下载及预处理

下载数据并转化为fastq格式：

#测试数据sra号
$ cat ERR.list
ERR1761568
ERR1761567
ERR1761566
ERR1761565
ERR1761564
ERR1761563
ERR1761562
ERR1761561
$ cat ERR.list | while read a;do echo "wget -c -nd -r -np -k -L -p -nd ftp://ftp.ncbi.nlm.nih.gov/sra/sra-instant/reads/ByRun/sra/ERR/${a:0:6}/$a/${a}.sra";done >run.down.sh
$ cat ERR.list | while read a;do echo "fastq-dump --split-3 ${a}.sra";done >run.fastqdump.sh

数据过滤：

$ cat ERR.list | while read a;do echo "trimmomatic PE -threads 4 ${a}_1.fastq ${a}_2.fastq ${a}_R1.PE.fastq ${a}_R1.SE.fastq ${a}_R2.PE.fastq ${a}_R2.SE.fastq SLIDINGWINDOW:4:15  LEADING:3 TRAILING:3 MINLEN:36 &>${a}.log";done >run.trim.sh

在NCBI上下载多个候选参考基因组数据：

#仅下载部分基因组数据
$ cat genome.list
Klebsiella_pneumoniae_KCTC_2242 ftp://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/000/220/485/GCF_000220485.1_ASM22048v1
Klebsiella_pneumoniae_ATCC_BAA-2146     ftp://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/000/364/385/GCF_000364385.3_ASM36438v3
Klebsiella_pneumoniae_500_1420  ftp://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/000/406/765/GCF_000406765.2_ASM40676v2
Klebsiella_pneumoniae_UHKPC33   ftp://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/000/417/085/GCF_000417085.2_ASM41708v2
Klebsiella_pneumoniae_DMC1097   ftp://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/000/417/225/GCF_000417225.2_ASM41722v2
Klebsiella_pneumoniae_UHKPC07   ftp://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/000/417/265/GCF_000417265.2_ASM41726v2
Klebsiella_pneumoniae_JM45              ftp://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/000/445/405/GCF_000445405.1_ASM44540v1
Klebsiella_pneumoniae_KP-1      ftp://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/000/465/975/GCF_000465975.2_ASM46597v2
Klebsiella_pneumoniae_CG43      ftp://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/000/474/015/GCF_000474015.1_ASM47401v1
Klebsiella_pneumoniae_PMK1      ftp://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/000/764/615/GCF_000764615.1_ASM76461v1
$ cat genome.list | while read a b;do echo "wget -c -r -nd -np -k -L -A fna.gz -P ${a} $b &>>down.log";done >run.down.sh

选择候选的参考基因组时需注意尽量选择组装完整的基因组，同时获得数据以后从基因组中去除质粒序列，例如：

$ fastaDeal.pl --pat  NC_022082.1 Klebsiella_pneumoniae_JM45.fna >tmp && mv tmp Klebsiella_pneumoniae_JM45.fna

同时下载了ERR1761507作为后续建树的外群，对于该数据的分析同上。

2. 参考基因组选择

$ refrank.py -f *PE.fastq -r *fna -t 10 &>refrank.log &
#只是为了跑个pipeline，所以只用了一个数据做测试
$ refrank.py -f ERR176156[78]* -r *fna -t 10 &>refrank.log &

得到一个refRanking.txt文件，这个文件里面有软件找出的参考基因组信息。

我仅用了部分SRA以及部分后续参考基因组，软件跑出来的结果是Klebsiella pneumoniae UHKPC33，就用该基因组作为后续分析的参考基因组了。

3. call SNP

使用文章中的流程call snp (Varscan + snpfilter.py)。

#以ERR1761561样本为例：
#建索引，生成mpileup文件
$ bwa index -p ref Klebsiella_pneumoniae_UHKPC33.fna &>/dev/null &
$ bwa mem -t 6 ref ERR1761561_R1.PE.fastq ERR1761561_R2.PE.fastq 2>/dev/null | samtools view -b -@ 4 -F 4 - | samtools sort -@ 4 - | samtools mpileup -q 1 -f Klebsiella_pneumoniae_UHKPC33.fna - >ERR1761561.mpileup 2>ERR1761561.log &
#VarScan软件call snp
$ java -jar /public/software/varscan/VarScan.v2.4.3.jar mpileup2snp ERR1761561.mpileup --output-vcf 1 >ERR1761561.snp.vcf 2>ERR1761561.log &
#基于参考序列生成含有snp的序列文件
$ perl chg.pl Klebsiella_pneumoniae_UHKPC33.fna ERR1761561.snp.vcf >Klebsiella_pneumoniae_UHKPC33.snp.fna
#过滤snp
$ snpfilter.py Klebsiella_pneumoniae_UHKPC33.fna Klebsiella_pneumoniae_UHKPC33.snp.fna

# 批处理
$ cat ERR.list | while read a;do echo "bwa mem -t 6 ref ${a}_R1.PE.fastq ${a}_R2.PE.fastq 2>/dev/null | samtools view -b -@ 4 -F 4 - | samtools sort -@ 4 - | samtools mpileup -q 1 -f Klebsiella_pneumoniae_UHKPC33.fna - >${a}.mpileup 2>/dev/null && java -jar /public/software/varscan/VarScan.v2.4.3.jar mpileup2snp ${a}.mpileup --output-vcf 1>${a}.snp.vcf 2>/dev/null";done >run.sh
$ cat ERR.list | while read a;do perl chg.pl Klebsiella_pneumoniae_UHKPC33.fna ${a}.snp.vcf >>all.snp.fa;done
#过滤snp，去gap、ambiguous base以及基因组中的重复区域和前噬菌体序列等。做这步分析之前需要使用phaster在线工具注释基因组中的噬菌体序列，并使用ICEberg下的ICEfinder在线工具注释细菌基因组中的整合子结合元件等
$ snpfilter.py Klebsiella_pneumoniae_UHKPC33.fna all.snp.fa --mask region.txt

使用phaster在线工具过滤噬菌体序列。phaster在线注释出参考基因组中的噬菌体区域：

460061  490463
1724712 1733196
2918297 2978044
3294633 3340680
3480582 3519282
3955718 4006647
4213795 4235686
4696467 4709499

使用ICEfinder在线工具注释细菌基因组的整合子结合元件（integrative and conjugative elements）：

#   Name    Location    Length/bp   Type
1   Region1 740482..796624  56143   Putative ICE with T4SS
2   Region2 1792161..1819646    27486   Putative IME

最后将phaster和ICEfinder的结果手动合并成一个文本文件，这里命名为region.txt。该文件包括两列，第一列为起始位点，第二列为终止位点，中间以制表符分割。

460061  490463
740482  796624
1724712 1733196
1792161 1819646
2918297 2978044
3294633 3340680
3480582 3519282
3955718 4006647
4213795 4235686
4696467 4709499

snpfilter.py处理结果会生成fa文件，即过滤后并串联后的序列，可以用于下游分析建树。将生成的snp串联后的fa文件重命名为SNPFILTER.fa
吐槽一下，VarScan和snpfilter.py两个软件的说明文档写的真心差，有些地方靠猜。vcf格式需要转化成带有snp的序列文件这一点文章和说明文档里面都没有，完全是多次试验自己猜出来的。

chg.pl自己写的，脚本如下：

#!/usr/bin/perl -w
use strict;

open FA,"$ARGV[0]";
open SNP,"$ARGV[1]";

my $file = "$ARGV[1]";
$file =~ s/\.snp\.vcf//;

my $seq;
while(<FA>){
        unless(/^>/){
                chomp;
                $seq .= $_;
        }
}

while(<SNP>){
        unless(/^#/){
                chomp;
                my @F = split /\t/,$_;
                my $snp = $F[4];
                my $locate = $F[1];
                substr($seq,$F[1]-1,1) = $snp;
        }
}
print ">$file\n$seq\n";

这个方法call snp不算完美，更棒的教程可参考Cosmika Goswami和Umer Zeeshan Ijaz 的call snp流程**
流程具体链接见“http://userweb.eng.gla.ac.uk/cosmika.goswami/snp_calling/SNPCalling.html”

4.构建系统发育树

#可以使用一下流程。但是测试的时候用的服务器内存有限，muscle这一步太吃内存了，改用了耗内存更小的mafft。
$ muscle -in SNPFILTER.fa -out SNPFILTER.aln && \
BeforePhylo.pl -type=dna -trim -conc=raxml -output=phylip SNPFILTER.aln && \
mv output.phy SNPFILTER.phy && \
raxmlHPC -T 6 -f a -x 12345 -p 12345 -# 1000 -m GAMMA -s SNPFILTER.phy -n SNPFILTER >/dev/null &

#改用一下流程
$ mafft --auto --phylipout --thread 10 --reorder SNPFILTER.fa > SNPFILTER.phy

BeforePhylo.pl这个脚本可以在这个路径下载：https://github.com/qiyunzhu/BeforePhylo

细菌全基因组测序数据构建SNP系统发育树

1. 数据下载及预处理

2. 参考基因组选择

3. call SNP

4.构建系统发育树

猜你喜欢

热点阅读