DADA2+Kraken2分析16s之细菌

DADA2+Kraken2分析16s之细菌

  DADA2在QIIME2中是以命令行形式呈现，有时候我们不想装QIIME2，但又想用命令行来运行DADA2。那可咋搞？我们可自己封装R代码为命令行，这个方式对于不熟悉代码的人来说效率比较低，我们可以从QIIME2里面把对应的DADA2运行的代码抠出来，点击我即可链接到地址；我们可以看到有run_dada_paired.R和run_dada_single.R两个文件，根据自己需要下载即可，这里以run_dada_paired.R为例。下载后，我们需要修改命令行传参的参数以及增加一个参数用来指定序列格式：

inp.dirF <- args[[1]]
inp.dirR <- args[[2]]
out.path <- args[[3]]
out.track <- args[[4]]
filtered.dirF <- args[[5]]
filtered.dirR <- args[[6]]
truncLenF <- as.integer(args[[7]])
truncLenR <- as.integer(args[[8]])
trimLeftF <- as.integer(args[[9]])
trimLeftR <- as.integer(args[[10]])
maxEEF <- as.numeric(args[[11]])
maxEER <- as.numeric(args[[12]])
truncQ <- as.integer(args[[13]])
poolMethod <- args[[14]]
chimeraMethod <- args[[15]]
minParentFold <- as.numeric(args[[16]])
nthreads <- as.integer(args[[17]])
nreads.learn <- as.integer(args[[18]])

以上内容可以更改为:

suppressMessages(library(optparse))


if (TRUE){
  option_list = list(
  
    make_option(c("-s", "--filetype"), type="character", default=".fq",
                help="filetype [default %default]"),
    make_option(c("-F", "--inp.dirF"), type="character", default=NULL,
                help="inp.dirF [default %default]"),
    make_option(c("-R", "--inp.dirR"), type="character", default=NULL,
                help="inp.dirR [default %default]"),
				
	make_option(c("-a", "--filtered.dirF"), type="character", default=NULL,
                help="filtered.dirF [default %default]"),
	make_option(c("-b", "--filtered.dirR"), type="character", default=NULL,
                help="filtered.dirR [default %default]"),				
				
	make_option(c("-c", "--truncLenF"), type="double", default=240,
                help="truncLenF [default %default]"),				
				
	make_option(c("-d", "--truncLenR"), type="double", default=160,
                help="truncLenR [default %default]"),			
				
				
	make_option(c("-e", "--trimLeftF "), type="double", default=0,
                help="trimLeftF  [default %default]"),			
							
	make_option(c("-g", "--trimLeftR"), type="double", default=0,
                help="trimLeftR [default %default]"),								
	make_option(c("-n", "--nthreads"), type="double", default=0,
                help="nthreads [default %default]"),			
										
	make_option(c("-u", "--out.track"), type="character", default=NULL,
                help="out dir <dir> [default %default]"),
	
    make_option(c("-o", "--out.path"), type="character", default=NULL,
                help="out dir <dir> [default %default]")
  )
  opts = parse_args(OptionParser(option_list=option_list))
}

inp.dirF <- opts$inp.dirF
inp.dirR <- opts$inp.dirR
out.path <- opts$out.path
out.track <- opts$out.track
filtered.dirF <- opts$filtered.dirF
filtered.dirR <- opts$filtered.dirR
truncLenF <- as.integer(opts$truncLenF)
truncLenR <- as.integer(opts$truncLenR)
trimLeftF <- as.integer(opts$trimLeftF)
trimLeftR <- as.integer(opts$trimLeftR)
out.track <- opts$out.track
##以下参数我设置为默认值了，可根据自己需要来增加命令行参数
maxEEF <- as.numeric(2)
maxEER <- as.numeric(2)
truncQ <- as.integer(2)
poolMethod <- "independent"
chimeraMethod <- "consensus"
minParentFold <- as.numeric(1)
nthreads <- as.integer(0)
nreads.learn <- as.integer(1000000)

上面给出了:

 make_option(c("-s", "--filetype"), type="character", default=".fq",
                help="filetype [default %default]"),

  这是用来指定序列的后缀，假如说样本名规则为：样本名.1.fq，那么-s参数对应的就输入：.fq，需要修改源码为:

unfiltsF <- list.files(inp.dirF, pattern=".fastq.gz$", full.names=TRUE)
unfiltsR <- list.files(inp.dirR, pattern=".fastq.gz$", full.names=TRUE)
filtsF <- list.files(filtered.dirF, pattern=".fastq.gz$", full.names=TRUE)
filtsR <- list.files(filtered.dirR, pattern=".fastq.gz$", full.names=TRUE)
##将pattern的内容修改为：paste(opts$filetype,"$",sep="")
unfiltsF <- list.files(inp.dirF, pattern=paste(opts$filetype,"$",sep=""), full.names=TRUE)
unfiltsR <- list.files(inp.dirR, pattern=paste(opts$filetype,"$",sep=""), full.names=TRUE)
filtsF <- list.files(filtered.dirF, pattern=paste(opts$filetype,"$",sep=""), full.names=TRUE)
filtsR <- list.files(filtered.dirR, pattern=paste(opts$filetype,"$",sep=""), full.names=TRUE)

  综上，即可运用命令行模式运行DADA2:

mrdir=/mnt/d/amnew
##${mrdir}/Rdada2为结果文件存放目录
mkdir  ${mrdir}/Rdada2
## ${mrdir}/temp/F和 ${mrdir}/temp/R分别存放过滤后的文件
mkdir -p  ${mrdir}/temp/F
mkdir -p  ${mrdir}/temp/R
##Run
time Rscript run_dada_paired.R -s ".fq" -F ${mrdir}/rawData/F -R ${mrdir}/rawData/R \
-a /mnt/d/amnew/temp/F -b ${mrdir}/temp/R  -e 25 -g 26 -c 260 -d 220 -n 0 \
-o ${mrdir}/Rdada2/output.tsv -u ${mrdir}/Rdada2/track.tsv

  output.tsv为每一ASV丰度表，第一列为序列内容，也是ASV的id，只需要把第一列提取出来，增加ASV的ID即可，output.tsv的结构如下:

#OTU ID	SC.1.fq	SCT.1.fq
GTGGGGAATATTGGACAATGGGCGCAAGCCTGATCCAGCCATGCCGCGTGGGTGAAGAAGGCCCTAGGGTTGTAAACCCCTTTCAGCGGGGAAGATAATGACGGTACCCGCAGAAGAAGCCCCGGCTAACTCCGTGCCAGCAGCCGCGGTAATACGGAGGGGGCTAGCGTTGTTCGGAATTACTGGGCGTAAAGCGCGCGTAGGCGGACCTGCAAGTCGGGGGTGAAAGCCCGAGGCTCAACCTCGGAACTGCCTTCGATACTGCGGGTCTCGAGTCCGGGAGAGGTGAGTGGAATTCCTAGTGTAGAGGTGAAATTCGTAGATATTAGGAAGAACACCAGTGGCGAAGGCGGCTCACTGGCCCGGTACTGACGCTGAGGCGCGAAAGCGTGGGGAGCAAACAGG	308	525
GTGGGGAATATTGGACAATGGGCGCAAGCCTGATCCAGCCATGCCGCGTGGGTGATGAAGGCCCTAGGGTTGTAAAGCCCTTTCAGCGGGGAAGATAATGACGGTACCCGCAGAAGAAGCCCCGGCTAACTTCGTGCCAGCAGCCGCGGTAATACGAAGGGGGCTAGCGTTGCTCGGAATTACTGGGCGTAAAGCGCACGTAGGCGGATTGTTAAGTCGGGGGTGAAATCCTGGAGCTCAACTCCAGAACTGCCTTCGATACTGGCGATCTCGAGTCCGGGAGAGGTGAGTGGAACTGCGAGTGTAGAGGTGAAATTCGTAGATATTCGCAAGAACACCAGTGGCGAAGGCGGCTCACTGGCCCGGTACTGACGCTGAGGTGCGAAAGCGTGGGGAGCAAACAGG	142	251	

  在这里我们把此表的第一列改为“ASV_当前行数”，输出为asvtable.txt，并且把第一行的“.1.fq”删掉：

sed 's/#OTU ID/ASV_ID/g' ${mrdir}/Rdada2/output.tsv | awk 'NR>1{OFS="\t";$1="ASV_"(NR-1);$0=$0}1' \
| sed '1 s/.1.fq//g' > ${mrdir}/Rdada2/asvtab.txt

  得到的asvtable.tx结构如下:

ASV_ID	SC	SCT	
ASV_1	308	525	
ASV_2	142	251
ASV_3	0	0		
ASV_4	144	215	
ASV_5	0	43	

  接下来是代表序列的整理，只需要把output.tsv的第一列提取出来，并且加上“>ASV_行号”即可，这样代表序列的ID与asvtable.txt的ID就一样了。

awk '{print $1}' ${mrdir}/Rdada2/output.tsv | sed '1d' | awk '{$0=">ASV_"NR"\n"$0}1' > ${mrdir}/Rdada2/rep-seqs.fasta

  rep-seqs.fasta的结构如下:

>ASV_1
GTGGGGAATATTGGACAATGGGCGCAAGCCTGATCCAGCCATGCCGCGTGGGTGAAGAAGGCCCTAGGGTTGTAAACCCCTTTCAGCGGGGAAGATAATGACGGTACCCGCAGAAGAAGCCCCGGCTAACTCCGTGCCAGCAGCCGCGGTAATACGGAGGGGGCTAGCGTTGTTCGGAATTACTGGGCGTAAAGCGCGCGTAGGCGGACCTGCAAGTCGGGGGTGAAAGCCCGAGGCTCAACCTCGGAACTGCCTTCGATACTGCGGGTCTCGAGTCCGGGAGAGGTGAGTGGAATTCCTAGTGTAGAGGTGAAATTCGTAGATATTAGGAAGAACACCAGTGGCGAAGGCGGCTCACTGGCCCGGTACTGACGCTGAGGCGCGAAAGCGTGGGGAGCAAACAGG
>ASV_2
GTGGGGAATATTGGACAATGGGCGCAAGCCTGATCCAGCCATGCCGCGTGGGTGATGAAGGCCCTAGGGTTGTAAAGCCCTTTCAGCGGGGAAGATAATGACGGTACCCGCAGAAGAAGCCCCGGCTAACTTCGTGCCAGCAGCCGCGGTAATACGAAGGGGGCTAGCGTTGCTCGGAATTACTGGGCGTAAAGCGCACGTAGGCGGATTGTTAAGTCGGGGGTGAAATCCTGGAGCTCAACTCCAGAACTGCCTTCGATACTGGCGATCTCGAGTCCGGGAGAGGTGAGTGGAACTGCGAGTGTAGAGGTGAAATTCGTAGATATTCGCAAGAACACCAGTGGCGAAGGCGGCTCACTGGCCCGGTACTGACGCTGAGGTGCGAAAGCGTGGGGAGCAAACAGG
>ASV_3
GTGGGGGATCTTGCGCTAATGCGCGAAAGCGTGACGCAGCGACGCCGCGTGGGGGATGAAGGCCCTCGGGTTGTAAACCCCTTTCGGCAGGGACGAAGGGAAACTGACGGTACCTGCAGAAGAAGTCCCGGCTAACTACGTGCCAGCAGCCGCGGTAACACGTAGGGGACGAGCGTTGTCCGGATTCATTGGGCGTAAAGCGCGCGTAGGCGGCTCGGAAAGTCGGGCGTGAAATGCCGGGGCTCAACCCCGGGACTGCGTCCGATACTTCCGAGCTGGAGGTAGGTAGGGGAGATCGGAATTCCTGGTGTAGCGGTGAAATGCGCAGATATCAGGAGGAACACCGGTGGCGAAGGCGGATCTCTGGGCCTTACCTGACGCTGAGGCGCGAAAGCTAGGGGAGCGAACAGG

  如上，就得到了ASV丰度表和代表序列，经测试，上述结果与在qiime2里跑的结果一致。
接下来用kraken2进行物种注释:

kraken2 --db  /mnt/d/canopy/silva138 --threads 16 \
--report ${mrdir}/Rdada2/SAMPLE.kreport2 ${mrdir}/Rdada2/rep-seqs.fasta > ${mrdir}/Rdada2/SAMPLE.kraken2

  SAMPLE.kraken2文件里包含了每一ASV对应的物种信息，接下来是ASV表抽平、物种丰度表定量，可移步到微信公众号去查看，点击我可跳转。有时候我们会删除那些在多个样本中为0的ASV，这时候用R或其他工具可快速完成。
  最后，我们创建了学术交流QQ群：335774366。欢迎有兴趣的朋友加入→指导。

说明，DADA2的代码来源于github上QIIME2项目的插件:

https://github.com/qiime2/q2-dada2/blob/master/q2_dada2/assets/run_dada_paired.R