由于沒有服務(wù)器力喷，只能自己搗鼓著用windows的cmd命令行來批量下載一些數(shù)據(jù)足删，但是windows下高通量測(cè)序數(shù)據(jù)（如sra）下載體驗(yàn)真的很差，沒辦法匣距，只能硬著頭皮去搞
一、下載sratoolkit軟件
一）軟件位置：https://trace.ncbi.nlm.nih.gov/Traces/sra/sra.cgi?view=software
圖1.1.1 選擇箭頭所指windows版本
二）軟件安裝
- 將下載好的sratoolkit軟件解壓縮哎壳，然后配置環(huán)境變量
圖1.2.1
- 打開windows環(huán)境變量設(shè)置窗口：按快捷鍵
win+R后,輸入sysdm.cpl,然后回車,完事毅待，具體步驟請(qǐng)跟著后續(xù)圖示來。
圖1.2.2
圖1.2.3
圖1.2.4
圖1.2.5
三）檢測(cè)軟件是否可以正常工作
-
快捷鍵Win+R归榕，輸入cmd進(jìn)入cmd命令行
圖1.3.1 進(jìn)入cmd命令行界面尸红，進(jìn)入軟件所在位置，如：
bio-tools\sratoolkit\sratoolkit.2.10.7-win64\bin\prefetch.exe -h刹泄，如成功顯示幫助信息外里，說明軟件安裝成功
圖1.3.2
二、下載測(cè)試數(shù)據(jù)
如果你自己沒有數(shù)據(jù)（
16S rRNA基因測(cè)序數(shù)據(jù)）特石，那就去下載別人上傳到網(wǎng)絡(luò)各大數(shù)據(jù)庫的數(shù)據(jù)盅蝗，比如以NCBI為例
-
打開NCBI官網(wǎng)：https://www.ncbi.nlm.nih.gov/，隨便輸入關(guān)鍵詞如
16S rRNA sequencing姆蘸，選擇SRA數(shù)據(jù)庫墩莫，進(jìn)入。
圖2.1.1 -
選擇所需下載的個(gè)體
圖2.2.1
圖2.2.2
- 獲得所下載個(gè)體的SRR號(hào)逞敷，然后用
sratoolkit軟件的prefetch函數(shù)批量下載狂秦，如下圖所示。
F:\>bio-tools\sratoolkit\sratoolkit.2.10.7-win64\bin\prefetch.exe --option-file C:\Users\jnzd_\Desktop\SRR_Acc_List.txt
可以看到下載信息：表示下載成功
2020-06-28T00:49:08 prefetch.2.10.7: 1) Downloading 'SRR12073641'...
2020-06-28T00:49:08 prefetch.2.10.7: Downloading via HTTPS...
2020-06-28T00:54:27 prefetch.2.10.7: HTTPS download succeed
2020-06-28T00:54:27 prefetch.2.10.7: verifying 'SRR12073641'...
2020-06-28T00:54:27 prefetch.2.10.7: 'SRR12073641' is valid
2020-06-28T00:54:27 prefetch.2.10.7: 1) 'SRR12073641' was downloaded successfully
2020-06-28T00:54:28 prefetch.2.10.7: 'SRR12073641' has 0 unresolved dependencies
2020-06-28T00:54:30 prefetch.2.10.7: 2) Downloading 'SRR12073642'...
2020-06-28T00:54:30 prefetch.2.10.7: Downloading via HTTPS...
2020-06-28T01:00:17 prefetch.2.10.7: HTTPS download succeed
2020-06-28T01:00:17 prefetch.2.10.7: verifying 'SRR12073642'...
2020-06-28T01:00:17 prefetch.2.10.7: 'SRR12073642' is valid
2020-06-28T01:00:17 prefetch.2.10.7: 2) 'SRR12073642' was downloaded successfully
2020-06-28T01:00:17 prefetch.2.10.7: 'SRR12073642' has 0 unresolved dependencies
2020-06-28T01:00:20 prefetch.2.10.7: 3) Downloading 'SRR12073643'...
2020-06-28T01:00:20 prefetch.2.10.7: Downloading via HTTPS...
2020-06-28T01:08:38 prefetch.2.10.7: HTTPS download succeed
2020-06-28T01:08:38 prefetch.2.10.7: verifying 'SRR12073643'...
2020-06-28T01:08:38 prefetch.2.10.7: 'SRR12073643' is valid
2020-06-28T01:08:38 prefetch.2.10.7: 3) 'SRR12073643' was downloaded successfully
2020-06-28T01:08:38 prefetch.2.10.7: 'SRR12073643' has 0 unresolved dependencies
2020-06-28T01:08:40 prefetch.2.10.7: 4) Downloading 'SRR12073644'...
2020-06-28T01:08:40 prefetch.2.10.7: Downloading via HTTPS...
2020-06-28T01:15:32 prefetch.2.10.7: HTTPS download succeed
2020-06-28T01:15:32 prefetch.2.10.7: verifying 'SRR12073644'...
2020-06-28T01:15:32 prefetch.2.10.7: 'SRR12073644' is valid
2020-06-28T01:15:32 prefetch.2.10.7: 4) 'SRR12073644' was downloaded successfully
2020-06-28T01:15:32 prefetch.2.10.7: 'SRR12073644' has 0 unresolved dependencies
- sra數(shù)據(jù)轉(zhuǎn)換為fastq格式
for /l %i in (41,1,44) do (fastq-dump.exe --split-files --gzip D:\R\R-exercise\ncbi-data\SRR120736%i\SRR120736%i.sra -O D:\R\R-exercise\ncbi-data\output)
# 本想用--split-3參數(shù)兰粉，但windows下一直報(bào)錯(cuò)故痊，不知道怎么回事，用s--split-files不會(huì)玖姑，-O接轉(zhuǎn)換后的數(shù)據(jù)存放目錄
圖2.2.3 轉(zhuǎn)換后的數(shù)據(jù)愕秫，以fastq結(jié)尾
三、R語言運(yùn)行data2
- 將轉(zhuǎn)換為fastq的16S rRNA 基因測(cè)序數(shù)據(jù)用data2進(jìn)行測(cè)序數(shù)據(jù)分析焰络，包括去除引物戴甩、去除低質(zhì)量數(shù)據(jù)、去重復(fù)闪彼、去嵌合體等步驟：
#############################################
## FGFP sample inference from amplicon data
## by Raul Tito & Rodrigo Bacigalupe
## date: March 28th 2018
#############################################
### Pipeline created from https://benjjneb.github.io/dada2/tutorial.html
### More information available at https://benjjneb.github.io/dada2/index.html
### This pipeline requires demultiplexed files: Two fastq files per sample (*R1 and *R2), they can be *.gz
### Load required modules
# module load R/3.4.2 ### very important to use this version
### Load necessary libraries and set seed
# 國(guó)內(nèi)用戶清華鏡像站甜孤，加速國(guó)內(nèi)用戶下載
site="https://mirrors.tuna.tsinghua.edu.cn/CRAN"
# 檢查是否存在Biocondoctor安裝工具，沒有則安裝
if (!requireNamespace("BiocManager", quietly = TRUE))
install.packages("BiocManager",repo=site)
# 加載安裝工具畏腕，并安裝dada2
library(BiocManager)
BiocManager::install("dada2")
library(dada2)
packageVersion("dada2")
set.seed(12345)
### Every new MiSeq or Hiseq run is quality checked for number of reads and presence of unused barcodes to estimate
### levels of carryover from prevous runs or contamination events. As part of this initial analysis demultipled files
### per each sample are generated (Using LotuS).
### Set directory and files
#' datasets used were downloaded from:https://trace.ncbi.nlm.nih.gov/Traces/study/?acc=SRP268421&o=acc_s%3Aa
path <- "D:\\R\\R-exercise\\ncbi-data\\output" # CHANGE to the directory containing your demultiplexed R1 and R2 fastq files
fns <- list.files(path) #列出某文件夾下的所有文件
#fns
####################
### load your data
####################
fastqs <- fns[grepl(".fastq.gz$", fns)] # grepl查找匹配后返回true或false缴川，而grep查找匹配后返回的是索引
fastqs <- sort(fastqs) # Sort ensures forward/reverse reads are in same order
### make sure that R1 is for forward read and R2 for reverse
fnFs <- fastqs[grepl(".1.fastq.gz", fastqs)] ## Just the forward read files
fnRs <- fastqs[grepl(".2.fastq.gz", fastqs)] ## Just the reverse read files
## Get sample names from the first part of the forward read filenames
sample.names <- sapply(strsplit(fnFs, ".1.fastq.gz"), `[`, 1) ## check if it is 1 or 2
## Fully specify the path for the fnFs and fnRs
fnFs <- file.path(path, fnFs)
fnRs <- file.path(path, fnRs)
- 雙端測(cè)序數(shù)據(jù)的質(zhì)量值作圖
###########################################
## Examine quality profiles of F & R reads
###########################################
pdf("plotQualityProfile.pdf", onefile=T)
plotQualityProfile(fnFs[1:2]) ## remove 20 plus 10 (primers and first odd bases)
plotQualityProfile(fnRs[1:2]) ## remove 20 plus 10 (primers and first odd bases)
dev.off()
圖3.1
- 根據(jù)質(zhì)量圖對(duì)數(shù)據(jù)進(jìn)行過濾和截取
##################################
## Perform filtering and trimming
##################################
filt_path <- file.path(path, "filtered") # Place filtered files in filtered/subdirectory
filtFs <- file.path(filt_path, paste0(sample.names, "_F_filt.fastq.gz"))
filtRs <- file.path(filt_path, paste0(sample.names, "_R_filt.fastq.gz"))
- 過濾雙端序列，F(xiàn)ilter the forward and reverse reads:Important to remove primers and low quality regions（去除引物和低質(zhì)量的數(shù)據(jù)）
out <- filterAndTrim(fnFs, filtFs, fnRs, filtRs, truncLen=c(130,200), ## Different settings tried,these are good for current primer constructs for MiSeq and HiSeq描馅，truncLen參數(shù)設(shè)定序列長(zhǎng)度范圍把夸，不在該范圍則過濾
trimLeft=c(30, 30), # 每條reads開頭的30個(gè)堿基一般質(zhì)量值偏低，建議過濾铭污，或者根據(jù)實(shí)際情況設(shè)定參數(shù)范圍恋日。
maxN=0, maxEE=c(2,2), truncQ=11, rm.phix=TRUE,
compress=TRUE, multithread=TRUE) #
head(out)
## Examine quality profiles of filtered reads
pdf("plotQualityProfile.filt.pdf", onefile=T)
plotQualityProfile(filtFs[1:2])
plotQualityProfile(filtRs[1:2])
dev.off()
圖3.2 過濾后的序列質(zhì)量分布圖
- 錯(cuò)誤率
## Learn the Error Rates
## DADA2算法使用機(jī)器學(xué)習(xí)構(gòu)建參數(shù)誤差模型膀篮，認(rèn)為每個(gè)擴(kuò)增子測(cè)序樣品都具有不同的誤差比率。該learnErrors方法通過交替估計(jì)錯(cuò)誤率和對(duì)參考樣本序列學(xué)習(xí)錯(cuò)誤模型岂膳，直到學(xué)習(xí)模型同真實(shí)錯(cuò)誤率收斂于一致誓竿。與許多機(jī)器學(xué)習(xí)方法一樣，算法有一個(gè)初始假設(shè)：數(shù)據(jù)中的最大可能錯(cuò)誤率就是只有最豐富的序列是正確的谈截，其余都是錯(cuò)誤筷屡。
#########################
## Learn forward error rates
errF <- learnErrors(filtFs, nread=1e6, multithread=TRUE) ## variable but this is the minimum number of reads
## Learn reverse error rates
errR <- learnErrors(filtRs, nread=1e6, multithread=TRUE) ## variable but this is the minimum number of reads
## Sample inference and merger of paired-end reads
mergers <- vector("list", length(sample.names))
names(mergers) <- sample.names
## Plot estimated error as sanity check
pdf("plotErrors_F.pdf", onefile=T)
plotErrors(errF, nominalQ=TRUE)
dev.off()
pdf("plotErrors_R.pdf", onefile=T)
plotErrors(errR, nominalQ=TRUE)
dev.off()
圖3.3 錯(cuò)誤率模型圖
- 去重復(fù)序列
與usearch去冗余步驟類似，僅僅保留重復(fù)序列中的一條序列
#########################
## Perform dereplication
#########################
## Dereplicate the filtered fastq files
derepRs <- derepFastq(filtRs, verbose=TRUE)
derepFs <- derepFastq(filtFs, verbose=TRUE)
# Name the derep-class objects by the sample names
names(derepFs) <- sample.names
names(derepRs) <- sample.names
####################
## Sample Inference
####################
## Apply the core sequence-variant inference algorithm to the dereplicated data
## Infer the sequence variants in each sample
## 基于錯(cuò)誤模型進(jìn)一步質(zhì)控
dadaFs <- dada(derepFs, err=errF, multithread=TRUE)
dadaRs <- dada(derepRs, err=errR, multithread=TRUE)
## Inspect the dada-class object returned by dada
dadaFs[[1]]
# 473 sequence variants were inferred from 5288 input unique sequences 從5288個(gè)輸入獨(dú)特序列中推斷出473個(gè)真實(shí)序列
- 去噪后的正反向序列進(jìn)行合并
合并的條件是：默認(rèn)情況下傻盟，僅當(dāng)正向和反向序列重疊至少
12個(gè)堿基并且在重疊區(qū)域中彼此相同時(shí)才輸出合并序列
## Merge the denoised forward and reverse reads
mergers <- mergePairs(dadaFs, derepFs, dadaRs, derepRs, verbose=TRUE)
## Inspect the merger data.frame from the first sample
head(mergers[[1]])
mergers對(duì)象格式為R語言的數(shù)據(jù)框速蕊，數(shù)據(jù)框中包含序列及其豐度信息嫂丙，未能成功合并的序列被刪除娘赴。
- 構(gòu)建ASV表
開始構(gòu)建 amplicon sequence variant（ASV）表，跟啤；類似于我們傳統(tǒng)的OTU表一樣诽表，使用makeSequenceTable參數(shù)進(jìn)行構(gòu)建
############################
## Construct sequence table
############################
seqtab <- makeSequenceTable(mergers)
## Get dimensions
dim(seqtab)
# [1] 4 1237 可以看到我們是4個(gè)樣本，然后最終得到1237個(gè)ASV
## Inspect distribution of sequence lengths
table(nchar(getSequences(seqtab)))
- 去除嵌合體序列
Dada核心質(zhì)控方法糾正了一些錯(cuò)誤隅肥，
但嵌合體仍然存在竿奏。幸運(yùn)的是，去噪后序列的準(zhǔn)確性使得識(shí)別嵌合體比處理模糊OTU時(shí)更簡(jiǎn)單腥放。
## Remove chimeras
###################
## Remove chimeric sequences:
seqtab.nochim <- removeBimeraDenovo(seqtab, method="consensus", multithread=TRUE, verbose=TRUE)
dim(seqtab.nochim)
sum(seqtab.nochim)/sum(seqtab)
- 統(tǒng)計(jì)序列信息
根據(jù)上述幾個(gè)步驟的序列篩選泛啸，來mark一下原始序列到最終可以采用的時(shí)候其艱辛旅程吧
####################################
## Track reads through the pipeline
####################################
getN <- function(x) sum(getUniques(x))
track <- cbind(out, sapply(dadaFs, getN), sapply(mergers, getN), rowSums(seqtab), rowSums(seqtab.nochim))
# If processing a single sample, remove the sapply calls: e.g. replace sapply(dadaFs, getN) with getN(dadaFs)
colnames(track) <- c("input", "filtered", "denoised", "merged", "tabled", "nonchim")
rownames(track) <- sample.names
head(track)
圖3.4 序列篩選過程統(tǒng)計(jì)
- 物種注釋
###################
## Assign taxonomy
###################
#' 幾個(gè)常用的注釋數(shù)據(jù)庫：https://benjjneb.github.io/dada2/training.html
taxHS <- assignTaxonomy(seqtab.nochim, "data/rdp_train_set_16.fa.gz", multithread=TRUE) ## CHANGE to directory and pertinent database
taxHS_spe <- addSpecies(taxHS, "data/rdp_species_assignment_16.fa.gz") # 需要注釋種水平另外加載相應(yīng)數(shù)據(jù)庫
colnames(taxHS) <- c("Kingdom", "Phylum", "Class", "Order", "Family", "Genus")
unname(head(taxHS)) # Remove the names or dimnames attribute of an R object when presentation
unname(tail(taxHS))
#######################
## Write data to files
#######################
write.table(track, file = "track.tsv", quote=FALSE)
seqtab.nochim <- as.data.frame(seqtab.nochim) # 將其轉(zhuǎn)為data.frame
names(seqtab.nochim) <- paste("ASV",seq(1,dim(seqtab.nochim)[2],1), sep = "_")
# ASV標(biāo)簽名稱太長(zhǎng)，這里將其轉(zhuǎn)換下命名
write.table(seqtab.nochim, file = "sequence_table_SV.tsv", quote=FALSE)
write.table(taxHS, file = "taxa_SV.tsv", quote=FALSE)
taxHS_spe <- as.data.frame(taxHS_spe)
rownames(taxHS_spe) <- paste("ASV",seq(1,dim(seqtab.nochim)[2],1), sep = "_")
write.table(taxHS_spe, file = "taxa_SV_spe.txt", quote = FALSE)
圖3.5 ASV物種注釋結(jié)果
圖3.6 個(gè)體對(duì)應(yīng)的ASV豐度表
參考
- qiime2官網(wǎng)資料https://docs.qiime2.org/2020.6/
- DADA2中文教程v1.8https://blog.csdn.net/woodcorpse/article/details/86773312
