本次分析使用到的数据及数据处理流程参考文献“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
选择候选的参考基因组时需注意尽量选择组装完整的基因组，同时获得数据以后从基因组中去除质粒序列，例如：
$ count.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
