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ファイルをエクスポートして取得します.または、次のコマンドでエクスポートします.結果は同じです.
上R言語
結果をざっと見てみると、主に1、2つの属が2つの高レベルの分類に属している場合があります.例えば、梭菌属などは、まずコードで結合し、フォーマットして命名する必要があります.また、このレベルに達していない分類もあります.少し処理する必要があります.私のコードは以下の通りです.間違いがあるかもしれません.交流を歓迎します.R言語のレベルが高くないため、必ず最適化の余地があり、指摘も歓迎します.
これで!
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)
これで!