Qiime2 数据导入

Posted 2023-04-30

参考技术A 在qiime2中进行任何数据分析的第一步永远是将数据导入qiime2并储存为qiime对象(.qza)。qiime2 接受导入的数据类型很多，包括从刚下机到分析过程中产生的任何常用数据格式和类型，如果遇到下面没有提到的数据类型或格式，可以去 QIIME 2 Forum 寻找帮助。

标准的EMP单端测序文件应该包括两个fastq.gz：测序reads和barcode reads。这种格式下的的序列是混合的，例如：所有样品的raw data混合在一个.gz文件中。而文件中的records顺序是联系barcode和测序结果的桥梁，也是分离混合数据的关键。

将两个.gz文件放到一个文件夹如：emp-single-end-sequences中。

标准的EMP双端测序文件应该包含三个fastq.gz：forward sequence reads, reverse sequence reads 和 barcode reads。这种格式下的的序列也是混合的，例如：所有样品的raw data混合在 forward 和 reverse .gz文件中。而文件中的records顺序是联系barcode和测序结果的桥梁，也是分离混合数据的关键。

将三个.gz文件放到一个文件夹如：emp-paired-end-sequences中。

Casava 1.8单端测序结果是一个.zip文件夹，里面包含一系列的fastq.gz文件，每一个sample对应一个.gz文件。

下划线分割的各区域从左往右包括：

与单端测序结果相似，但对任一sample，双端测序结果包括两个fastq.gz文件，由R1和R2区分read 方向。

如果测序结果既不是EMP又不是Casava格式，我们就需要先自己创建一个 “manifest file”， 然后再用import 工具导入。

Fastq数据有 四种常用格式变体 ，导入时必须在--input-format 中指定。下面提供SingleEndFastqManifestPhred33V2 的导入方法，其他类似。

qiime2目前支持seqs.fna文件导入，该文件每个record都有两行：header 和 sequence的fasta文件；每条序列只能是一行，不能拆分为多行；每条序列的ID必须遵循 <sample-id>_<seq-id> 格式。 <sample-id> 是序列所属样本的标识符， <seq-id> 是其样本中序列的标识符。

这类数据结构上和上面的seqs.fna文件相似，但它是unaligned (i.e., do not contain - or . characters) ，并且包含未知序列(N)，有些插件不支持分析含N的序列。:

该类文件与上述的seqs.fna文件相似，但其序列是一一对齐并长度相同的，同样也可以包含未知碱基N，有些插件不支持分析含N的序列。

qiime2 支持 newick 树文件输入。

遇到的问题很可能以前有人已经遇到过，并且在获得了很好的技术支持并解决了。如果没有，上面也有很多热心的人提供帮助和思路。

使用R语言获得16S物种丰度

 还是获得16S物种丰度得老问题，最近在一台新机器上安装qiime1，发现有报错，对于这种停止维护的软件，也是正常现象吧，于是想别的办法解决，恰巧最近读R几本R语言的入门书，发现prop.table()这个函数是可以实现相关功能的，于是学习使用下。可能你早已会做这个啦，还是分享一下，看看有没有人需要。

从qiime2的结果开始

qiime2可以按级别导出物种的计数，通过view.qiime2.org或者view.qiime2.cn（后者是我自己镜像的前面一个网站） 查看taxa-bar-plots.qzv导出csv文件获得。或者通过下面的命令导出，结果应该是相同的。

# 按门和属水平合并，并统计
#门
qiime taxa collapse \
  --i-table table_final1.qza \
  --i-taxonomy taxonomy.qza \
  --p-level 2 \
  --o-collapsed-table table-l2.qza
#科
qiime taxa collapse \
  --i-table table_final1.qza \
  --i-taxonomy taxonomy.qza \
  --p-level 5 \
  --o-collapsed-table table-l5.qza
#属
qiime taxa collapse \
  --i-table table_final1.qza \
  --i-taxonomy taxonomy.qza \
  --p-level 6 \
  --o-collapsed-table table-l6.qza
#species
qiime taxa collapse \
  --i-table table_final1.qza \
  --i-taxonomy taxonomy.qza \
  --p-level 7 \
  --o-collapsed-table table-l7.qza
#导出 tsv
for file in ./table-l*.qza
do
  base=$(basename $file .qza)
  echo $base
  qiime tools export \
  --input-path $file \
  --output-path $base
  biom convert --to-tsv -i $base/feature-table.biom -o $base.tsv
done

上R语言

粗略看了下结果，主要是有一两个属竟然分属于两个高级别分类的情况，比如梭菌属等，需要先用代码合并下，并格式化命名，还有就是有的分类没有到达这个级别的，也需要稍作处理。我的代码如下，可能存在错误，欢迎交流指正。由于R语言水平不高，一定有优化的余地，也欢迎指出。

#读取文件
sample <- paste('table-l6','.tsv', sep = "")
df <- read.table(sample, header = TRUE, sep = '\t',comment.char="",skip=1, stringsAsFactors = FALSE)
#整理属名，去除多余信息
genus <- c()
for (j in 1:length(df[, 1])) {
  p <- paste(strsplit(df[, 1][j], "__")[[1]][7], j)
  if(strsplit(p, ' ')[[1]][1]=="NA"){
    if(strsplit(df[, 1][j], "__")[[1]][6]==';'){
      p <- df[, 1][j]
    }
    else if(strsplit(df[, 1][j], "__")[[1]][6]==';g'){
      p <- df[, 1][j]
    }
    else if(strsplit(p, ' ')[[1]][1]=="NA")  p <- paste(strsplit(strsplit(df[, 1][j], "__")[[1]][6], ';g')[[1]][1], j)
    else {p <-  paste(strsplit(strsplit(df[, 1][j], "__")[[1]][6], ';')[[1]][1], j) }
  }
  genus <- c(genus, p)
}
#属重命名
df <- df[-1]
row.names(df) <- genus
#合并相同属
get_genus_summary <- function(df, bacterium){
  Bac_name <- df[grepl(bacterium, rownames(df)),]
  Bac_name <- sapply(Bac_name, as.numeric)
  bac <- colSums(Bac_name)
}
df_new <- data.frame()
for (bact in 1:length(row.names(df))) {
  if(grepl("\\[",row.names(df)[bact])) {
    bac_name <- strsplit(strsplit(strsplit(row.names(df)[bact], ' ')[[1]][1], '\\[')[[1]][2], "\\]")[[1]][1]
  }else {
    bac_name <- strsplit(row.names(df)[bact], ' ')[[1]][1]
  }
 
  if (length(row.names(df[grepl(bac_name, rownames(df)),])) > 1) {
    if(!(bac_name %in% row.names(df_new))) {
      bac <- get_genus_summary(df, bac_name)
      df_new <- rbind(df_new, bac)
      row.names(df_new)[length(row.names(df_new))] <- bac_name
    }else next
  }
  else {
    df_new <- rbind(df_new, df[bact,])
    row.names(df_new)[length(row.names(df_new))] <- bac_name
  }
}
#获得比例数据，先转化成矩阵，2代表以列求和
df_new  <- prop.table(data.matrix(df_new), 2)

就这样啦！

- END -