-
Notifications
You must be signed in to change notification settings - Fork 0
/
Decontam_Batches_Tutorial.R
32 lines (26 loc) · 1.17 KB
/
Decontam_Batches_Tutorial.R
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
asv_df <- read.delim("/Users/jrabasc/Desktop/github/Decontam_Validation/dada2_output/dada2_feature_table/table.tsv", check.names=FALSE)
rownames(asv_df) <- asv_df[, 1]
asv_df <- asv_df[, -1]
asv_df<-t(asv_df)
numero_df <- as.matrix(sapply(asv_df, as.numeric))
metadata_df<-read.csv(file = "/Users/jrabasc/Desktop/github/Decontam_Validation/dada2_output/metadata_non_HMS_samples_removed.csv", check.names=FALSE)
rownames(metadata_df) <- metadata_df[, 1]
meta_data_cols <-function(asv_df, metadata_df, control.col){
control_vec<-c()
index<-0
for (id in colnames(metadata_df)) {
index=index+1
if(id == control.col){
control_vec<-metadata_df[,c(index)]
}
}
# We need to align the metadata to the asv table
mapped_ids <- match( rownames(asv_df),rownames(metadata_df))
control_vec <- control_vec[mapped_ids]
# drop sample IDs not in the asv table
control_vec <- na.omit(control_vec)
return(control_vec)
}
control_vec <- meta_data_cols(asv_df, metadata_df, "exp_sample_type")
true_false_control_vec<-grepl("control",control_vec)
prev_contam <- isContaminant(asv_df, neg=true_false_control_vec, threshold=0.1, detailed=TRUE, normalize=TRUE, method='prevalence')