R言語で16 S種の豊富さを得る

17079 ワード

やはり16 S種の豊富さを得るのが古い問題で、最近1台の新しい機械の上でqiime 1をインストールして、間違いを発見して、このようなメンテナンスを停止するソフトウェアに対して、正常な現象でしょう、そこで別の方法で解決して、ちょうど最近RのいくつかのR言語の入門書を読んで、propを発見します.table()という関数は関連機能を実現できるので,学習して使用する.あなたはとっくにこれをすることができるかもしれませんが、やはり分かち合って、誰かが必要かどうかを見てみましょう.
qiime 2の結果から
qiime 2は、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言語
結果をざっと見てみると、主に1、2つの属が2つの高レベルの分類に属している場合があります.例えば、梭菌属などは、まずコードで結合し、フォーマットして命名する必要があります.また、このレベルに達していない分類もあります.少し処理する必要があります.私のコードは以下の通りです.間違いがあるかもしれません.交流を歓迎します.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 p
    }
    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)

これで!