From epidome to treeHeatmap

Leave a reply

  1. (namely epidome-processing step 2.6) Classification using 98% similarity.

    #Input Files: yycH_seqtab_nochim.csv and ../epidome/DB/yycH_ref_aln.fasta
#Output Files: yycH_seqtab_ASV_seqs.fasta, yycH_seqtab_ASV_blast.txt, and yycH_seqtab.csv.classified.csv

python3 ../epidome/scripts/ASV_blast_classification.py yycH_seqtab_nochim.csv yycH_seqtab_ASV_seqs.fasta ../epidome/DB/yycH_ref_aln.fasta yycH_seqtab_ASV_blast.txt yycH_seqtab.csv.classified.csv 98
python3 ../epidome/scripts/ASV_blast_classification.py g216_seqtab_nochim.csv g216_seqtab_ASV_seqs.fasta ../epidome/DB/g216_ref_aln.fasta g216_seqtab_ASV_blast.txt g216_seqtab.csv.classified.csv 98

  2. Manually check the reults, ensuring no seq34;seq35 occurs.

    #Delete all records in *_seqtab.csv.classified.csv the alignment length < 448 (delete records 110 and 111) in g216, or < 440 (del record 90) in yycH.

#Confirm the results
jhuang@hamburg:/media/jhuang/Elements/Data_Luise_Epidome_batch3/run_2023_98$ cut -d$'\t' -f4 g216_seqtab_ASV_blast.txt | sort -u
446
447
448
jhuang@hamburg:/media/jhuang/Elements/Data_Luise_Epidome_batch3/run_2023_98$ cut -d$'\t' -f4 yycH_seqtab_ASV_blast.txt | sort -u
475
476

#It is not neccesary to delete correponding records 90 in yycH_* and g216_seqtab_ASV_blast.txt, since the two files are not used. 
#sed -i -e s/seq//g yycH_seqtab_ASV_blast.txt   #length=476
#sed -i -e s/seq//g g216_seqtab_ASV_blast.txt   #length=448
##;-->""
#sed -i -e s/';'//g yycH_seqtab_ASV_blast.txt
#sed -i -e s/';'//g g216_seqtab_ASV_blast.txt

cp g216_seqtab.csv.classified.csv g216_seqtab.csv.classified.backup
cp yycH_seqtab.csv.classified.csv yycH_seqtab.csv.classified.backup

# (Option 1): Processing *.classified.csv, ensuring the format seq24,21 or seq24.
#https://github.com/ssi-dk/epidome/blob/master/example_data/190920_run1_G216_seqtab_from_dada2.csv.classified.csv
##DEBUG using LibreOffice, e.g. libreoffice --calc yycH_seqtab.csv.classified_noNA.csv after adding "ID"; at the corner.
#Final goal: "seqseq31;,seq30;" --> "seq31,30"
#;,seq --> ,
sed -i -e s/";,seq"/","/g g216_seqtab.csv.classified.csv
sed -i -e s/";,seq"/","/g yycH_seqtab.csv.classified.csv
#;"; --> ";
sed -i -e s/";\";"/"\";"/g g216_seqtab.csv.classified.csv
sed -i -e s/";\";"/"\";"/g yycH_seqtab.csv.classified.csv
sed -i -e s/seqseq/seq/g g216_seqtab.csv.classified.csv
sed -i -e s/seqseq/seq/g yycH_seqtab.csv.classified.csv
#confirm!
cut -d';' -f3 g216_seqtab.csv.classified.csv
cut -d';' -f3 yycH_seqtab.csv.classified.csv

# (Option 2): To reduce the unclassified, rename seq31,seq30 --> seq31 in g216,  seq37,seq36 --> seq36 in yycH.
sed -i -e s/"seq40,seq31,seq30"/"seq40"/g g216_seqtab.csv.classified.csv
sed -i -e s/"seq31,seq30"/"seq31"/g g216_seqtab.csv.classified.csv
sed -i -e s/"seq5,seq2"/"seq2"/g        sed -i -e s/"seq40,seq31,seq30"/"seq40"/g g216_seqtab.csv.classified.csv
sed -i -e s/"seq31,seq30"/"seq31"/g g216_seqtab.csv.classified.csv
sed -i -e s/"seq5,seq2"/"seq2"/g g216_seqtab.csv.classified.csv
grep ",seq" g216_seqtab.csv.classified.csv  #should return no record.
#sed -i -e s/"seq22,20"/"seq20"/g g216_seqtab.csv.classified.csv
#sed -i -e s/"seq24,21"/"seq21"/g g216_seqtab.csv.classified.csv
#sed -i -e s/"seq21,17"/"seq21"/g g216_seqtab.csv.classified.csv
> epidome_object$p1_seqs    #vim g216_seqtab.csv.classified_noNA.csv
  [1] "seq28"       "seq5"        "seq40"       "seq31,30"(seq31)     "seq14"      
  [6] "seq20"       "seq22"       "seq26"       "seq37"       "seq21"      
[11] "seq29"       "seq24"       "seq1"        "seq26"       "seq40"      
[16] "seq21"       "seq18"       "seq40"       "seq16"       "seq37"      
[21] "seq5,2"      "seq18"       "seq11"       "seq14"       "seq26"      
[26] "seq31,30"(seq31)     "seq37"       "seq26"       "seq1"        "seq21"      
[31] "seq3"        "seq40"       "seq26"       "seq28"       "seq40"      
[36] "seq10"       "seq5"        "seq40,31,30"(seq40)  "seq37"       "seq37"      
[41] "seq37"       "seq40"       "seq28"       "seq26"       "seq37"      
[46] "seq22"       "seq28"       "seq19"       "seq5,2"      "seq26"      
[51] "seq40"       "seq40"       "seq26"       "seq40"       "seq28"      
[56] "seq28"       "seq28"       "seq37"       "seq26"       "seq40"      
[61] "seq26"       "seq37"       "seq28"       "seq28"       "seq1"       
[66] "seq28"       "seq5,2"(seq5)      "seq40"       "seq21"       "seq28"      
[71] "seq29"       
#"seq37"       "seq22"       "seq22,20"(seq20)    "seq31,30"(seq31)    
[76] "seq22"       "seq20"       "seq28"       "seq37"       "seq14"      
[81] "seq28"       "seq20"       "seq40"       "seq28"       "seq24,21" 
[86] "seq40"       "seq14"       "seq31,30"(seq31)      "seq1"        "seq1"       
[91] "seq31,30"(seq31)     "seq1"        "seq40"       "seq20"       "seq1"       
[96] "seq21"       "seq26"       "seq20"       "seq20"       "seq21,17"   
[101] "seq37"       "seq26"       "seq31,30"(seq31)      "seq31,30"(seq31)    "seq40"      
[106] "seq5"        "seq12"             

sed -i -e s/"seq37,seq36"/"seq36"/g yycH_seqtab.csv.classified.csv
sed -i -e s/"seq58,seq25,seq21,seq10"/"seq21"/g yycH_seqtab.csv.classified.csv
sed -i -e s/"seq23,seq22"/"seq23"/g yycH_seqtab.csv.classified.csv
sed -i -e s/"seq65,seq63,seq61,seq60"/"seq63"/g yycH_seqtab.csv.classified.csv
sed -i -e s/"seq65,seq63,seq53"/"seq63"/g yycH_seqtab.csv.classified.csv
sed -i -e s/"seq61,seq55"/"seq55"/g yycH_seqtab.csv.classified.csv
sed -i -e s/"seq34,seq28,seq25,seq10,seq9"/"seq34"/g yycH_seqtab.csv.classified.csv
sed -i -e s/"seq25,seq10,seq9"/"seq9"/g yycH_seqtab.csv.classified.csv
sed -i -e s/"seq25,seq10"/"seq10"/g yycH_seqtab.csv.classified.csv
grep ",seq" yycH_seqtab.csv.classified.csv  #should return no record.
> epidome_object$p2_seqs
[1] "34"            "43"            "seq37,seq36"-->36         "42"           
[5] "26"            "63"            "3"             "2"            
[9] "11"            "seq65,seq63,seq53"-->63      "14"            "38"           
[13] "55"            "3"             "65"            "seq58,seq25,seq21,seq10"-->21  
[17] "seq23,seq22"-->23         "53"            "56"            "34"           
[21] "62"            "14"            "35"            "seq61,seq55"-->55        
[25] "52"            "seq25,seq10"-->10         "seq65,seq63,seq61,seq60"-->63   "21"           
[29] "55"            "15"            "34"            "49"           
[33] "19"            "63"            "63"            "34"           
[37] "43"            "12"            "34"            "seq34,seq28,seq25,seq10,seq9"-->34
[41] "43"            "61"            "25"            "34"           
[45] "34"            "5"             "18"            "43"           
[49] "62"            "34"            "19"            "63"           
[53] "7"             
--------------
"seq37,seq36"-->36         "seq25,seq10,seq9"-->9       "56"           
[57] "11"            "49"            "10"            "63"           
[61] "seq37,seq36"-->36         "14"            "55"  
#-------------
"55"           
[65] "14"            "37,36"         "43"            "45"           
[69] "34"            "34"            "43"            "43"           
[73] "25"            "14"            "65,64,61,56"   "15,9"         
[77] "43"            "38"            "43"            "53"           
[81] "11"            "9"             "34"            "15"           
[85] "9"             "23,22"         "62"            "43"           
[89] "37,36"

# Common step after option 1 or option 2.
grep -v ";NA;" g216_seqtab.csv.classified.csv > g216_seqtab.csv.classified_noNA.csv
grep -v ";NA;" yycH_seqtab.csv.classified.csv > yycH_seqtab.csv.classified_noNA.csv

#filtering the small ASV in which the number < 500: INPUT: *_noNA.csv; OUTPUT: *_seqtab.csv.classified_filtered.csv
python3 g216_filtering.py  
    import pandas as pd
    # Load the file
    data = pd.read_csv('g216_seqtab.csv.classified_noNA.csv', delimiter=';')
    # Compute row sums excluding the first 2 columns
    data['row_sum'] = data.iloc[:, 2:].sum(axis=1)
    print(data.iloc[:, 2:])
    print(data.iloc[:, 2:].sum(axis=1))
    # Filter out rows where the sum is less than 100
    filtered_data = data[data['row_sum'] >= 500].drop(columns=['row_sum'])
    # Save the resulting dataframe to a new file
    filtered_data.to_csv('g216_seqtab.csv.classified_filtered.csv', sep=';', index=False)
python3 yycH_filtering.py

  3. draw plot from three amplicons: cutadapted_g216, cutadapted_yycH

    #Taxonomic database setup and classification
#- Custom databases of all unique g216 and yycH target sequences can be found at https://github.com/ssi-dk/epidome/tree/master/DB. 
#- We formatted our g216 and yycH gene databases to be compatible with DADA2’s assign-Taxonomy function and used it to classify the S. epidermidis ASVs with the RDP naive Bayesian classifier method (https://github.com/ssi-dk/epidome/tree/master/scripts).
#- ST classification of samples was performed using the g216 target sequence as the primary identifier. 
#- All g216 sequences unique to a single clonal cluster in the database were immediately classified as the matching clone, and in cases were the g216 sequence matched multiple clones, the secondary yycH target sequences were parsed to determine which clone was present. When this classification failed to resolve due to multiple potential combinations of sequences, ASVs were categorized as “Unclassified”. Similarly, g216 sequences not found in the database were labelled as “Novel”. 

#install.packages("pls")
#library(pls)
#install.packages("reshape")
#library(reshape)
#install.packages("vegan")
library('vegan') 
library(scales)
library(RColorBrewer)

#~/Tools/csv2xls-0.4/csv_to_xls.py g216_seqtab.csv.classified_noNA.csv yycH_seqtab.csv.classified_noNA.csv -d$'\t' -o counts.xls

setwd("/media/jhuang/Elements/Data_Luise_Epidome_batch3/run_2023_98_2")
#under r4-base
source("../epidome/scripts/epidome_functions.R")

ST_amplicon_table = read.table("../epidome/DB/epidome_ST_amplicon_frequencies.txt",sep = "\t")

epi01_table = read.table("g216_seqtab.csv.classified_filtered.csv",sep = ";",header=TRUE,row.names=1)
epi02_table = read.table("yycH_seqtab.csv.classified_filtered.csv",sep = ";",header=TRUE,row.names=1)
#> sum(epi01_table$AH.LH)
#[1] 78872
#> sum(epi02_table$AH.LH)
#[1] 86949

metadata_table = read.table("../metadata.txt",header=TRUE,row.names=1)
metadata_table$patient.ID <- factor(metadata_table$patient.ID, levels=c("P7","P8","P9","P10","P11","P12","P13","P14","P15","P16","P17","P18","P19","P20", "AH","AL","CB","HR","KK",  "MC","MR","PT2","RP","SA","XN"))

#epidome_object = setup_epidome_object(epi01_table,epi02_table, metadata_table=metadata_table)
primer1_table <- epi01_table
primer2_table <- epi02_table
#setup_epidome_object <- function(primer1_table,primer2_table,metadata_table) {
  #DEBUG: 3:ncol --> 2:ncol
  primer1_counts = primer1_table[,2:ncol(primer1_table)]
  primer2_counts = primer2_table[,2:ncol(primer2_table)]
  primer1_all_sample_names = colnames(primer1_counts)
  primer2_all_sample_names = colnames(primer2_counts)
  samples_with_both_primers = primer1_all_sample_names[which(primer1_all_sample_names %in% primer2_all_sample_names)]
  samples_missing_primer1_data = primer2_all_sample_names[which(!primer2_all_sample_names %in% primer1_all_sample_names)]
  samples_missing_primer2_data = primer1_all_sample_names[which(!primer1_all_sample_names %in% primer2_all_sample_names)]
  primer1_seqs = as.vector(primer1_table$Seq_number)
  primer2_seqs = as.vector(primer2_table$Seq_number)
  primer1_seqs[which(is.na(primer1_seqs))] = "seqUnclassified"
  primer2_seqs[which(is.na(primer2_seqs))] = "seqUnclassified"
  primer1_counts = primer1_table[,which(colnames(primer1_table) %in% samples_with_both_primers)]
  primer2_counts = primer2_table[,which(colnames(primer2_table) %in% samples_with_both_primers)]
  if (!missing(metadata_table)) {
    metadata_names = rownames(metadata_table)
    samples_missing_metadata = samples_with_both_primers[which(!samples_with_both_primers %in% metadata_names)]
    samples_missing_primer_data = metadata_names[which(!metadata_names %in% samples_with_both_primers)]
    include_names = metadata_names[which(metadata_names %in% samples_with_both_primers)]
    metadata_include = match(include_names,metadata_names)
    metadata_include = metadata_include[which(!is.na(metadata_include))]
    metadata_table = metadata_table[metadata_include,]
    #primer1_include = match(colnames(primer1_counts),include_names)
    primer1_include = match(include_names,colnames(primer1_counts))
    primer1_include = primer1_include[which(!is.na(primer1_include))]
    #primer2_include = match(colnames(primer2_counts),include_names)
    primer2_include = match(include_names,colnames(primer2_counts))
    primer2_include = primer2_include[which(!is.na(primer2_include))]
    epi1_table = primer1_counts[,primer1_include]
    epi2_table = primer2_counts[,primer2_include]
    epidome_object = list('p1_seqs'=primer1_seqs,'p1_table'=epi1_table,'p2_seqs'=primer2_seqs,'p2_table'=epi2_table,'metadata'=metadata_table,'sample_names'=include_names,'meta_variables'=colnames(metadata_table))
    print(paste0("Metadata loaded with ",length(metadata_names)," samples and ",ncol(metadata_table)," variables"))
    print(paste0(length(samples_missing_metadata)," samples are found in both primer sets but not in metadata: ",paste0(samples_missing_metadata,collapse = " ")))
    print(paste0(length(samples_missing_primer_data)," samples are found in metadata but is missing from one or both primer sets: ",paste0(samples_missing_primer_data,collapse = " ")))
    print(paste0(length(include_names)," samples are found in metadata and both tables and are included in epidome object"))

  } else {
    epi1_table = primer1_counts[,match(colnames(primer1_counts),samples_with_both_primers)]
    epi2_table = primer2_counts[,match(colnames(primer2_counts),samples_with_both_primers)]
    epidome_object = list('p1_seqs'=primer1_seqs,'p1_table'=epi1_table,'p2_seqs'=primer2_seqs,'p2_table'=epi2_table,'sample_names'=samples_with_both_primers)
    print(paste0("No metadata loaded"))
    print(paste0(length(samples_missing_primer2_data)," samples are found in p1 table but not in p2 table: ",paste0(samples_missing_primer2_data,collapse = " ")))
    print(paste0(length(samples_missing_primer1_data)," samples are found in p2 table but not in p1 table: ",paste0(samples_missing_primer1_data,collapse = " ")))
    print(paste0(length(samples_with_both_primers)," samples are found in both tables and are included in epidome object"))
  }
  #return(epidome_object)
#}

str(epidome_object)
#List of 7
# $ p1_seqs       : chr [1:57] "seq28" "seq5" "seq40" "seq31,30" ...
# $ p1_table      :'data.frame':    57 obs. of  51 variables:
# $ p2_seqs       : chr [1:53] "seq34" "seq43" "seq37,36" "seq42" ...
# $ p2_table      :'data.frame':    53 obs. of  51 variables:
# $ metadata      :'data.frame':    51 obs. of  4 variables:
# $ sample_names  : chr [1:51] "P7.Nose" "P8.Nose" "P9.Nose" "P10.Nose" ...
# $ meta_variables: chr [1:4] "patient.ID" "sample.site" "sample.type" "patient.sample.site"
unique(sort(epidome_object$p1_seqs))  # 57
#[1] "seq28"       "seq5"        "seq40"       "seq31,30"    "seq14"      
#[6] "seq20"       "seq22"       "seq26"       "seq37"       "seq21"      
#[11] "seq29"       "seq24"       "seq1"        "seq26"       "seq40"      
#[16] "seq21"       "seq18"       "seq40"       "seq16"       "seq37"      
#[21] "seq5,2"      "seq18"       "seq11"       "seq14"       "seq26"      
#[26] "seq31,30"    "seq37"       "seq26"       "seq1"        "seq21"      
#[31] "seq3"        "seq40"       "seq26"       "seq28"       "seq40"      
#[36] "seq10"       "seq5"        "seq40,31,30" "seq37"       "seq37"      
#[41] "seq37"       "seq40"       "seq28"       "seq26"       "seq37"      
#[46] "seq22"       "seq28"       "seq19"       "seq5,2"      "seq26"      
#[51] "seq40"       "seq40"       "seq26"       "seq40"       "seq28"      
#[56] "seq28"       "seq28"
unique(sort(epidome_object$p2_seqs))

#Image1
primer_compare = compare_primer_output(epidome_object,color_variable = "sample.type")
png("image1.png")
primer_compare$plot
dev.off()

eo_ASV_combined = combine_ASVs_epidome(epidome_object)
#> str(eo_ASV_combined)
#List of 7
#$ p1_seqs       : chr [1:21] "seq28" "seq5" "seq40" "seq31,30" ...
#$ p1_table      :'data.frame': 21 obs. of  51 variables:
#> eo_ASV_combined$p1_seqs  #21
#[1] "seq28"       "seq5"        "seq40"       "seq31,30"    "seq14"      
#[6] "seq20"       "seq22"       "seq26"       "seq37"       "seq21"      
#[11] "seq29"       "seq24"       "seq1"        "seq18"       "seq16"      
#[16] "seq5,2"      "seq11"       "seq3"        "seq10"       "seq40,31,30"
#[21] "seq19"
#eo_filtered = filter_lowcount_samples_epidome(eo_ASV_combined,500,500)

#Note that we run the following code instead run the method from the script due to modified code: 
count_table_checked = classify_epidome(eo_ASV_combined,ST_amplicon_table)
write.csv(count_table_checked, file="count_table_checked.csv")
strict_classifier=FALSE
#classify_epidome = function(epidome_object,ST_amplicon_table,strict_classifier=FALSE) {
#Adapt the code to the true epi01 and epi02 values, please give the complete code:
epidome_object_norm = normalize_epidome_object(eo_ASV_combined)
p1_seqs = unlist(lapply(as.vector(eo_ASV_combined$p1_seqs),function(x) substr(x,4,nchar(x))))
#> p1_seqs
#[1] "28"       "5"        "40"       "31,30"    "14"       "20"      
#[7] "22"       "26"       "37"       "21"       "29"       "24"      
#[13] "1"        "26"       "40"       "21"       "18"       "40"      
#[19] "16"       "37"       "5,2"      "18"       "11"       "14"      
#[25] "26"       "31,30"    "37"       "26"       "1"        "21"      
#[31] "3"        "40"       "26"       "28"       "40"       "10"      
#[37] "5"        "40,31,30" "37"       "37"       "37"       "40"      
#[43] "28"       "26"       "37"       "22"       "28"       "19"      
#[49] "5,2"      "26"       "40"       "40"       "26"       "40"      
#[55] "28"       "28"       "28" 
p2_seqs = unlist(lapply(as.vector(eo_ASV_combined$p2_seqs),function(x) substr(x,4,nchar(x))))
n_samples = length(eo_ASV_combined$sample_names)
n_p1_seqs = length(p1_seqs)
count_table = matrix(nrow = 0, ncol = n_samples,dimnames = list('ST'=c(),'Samples'=eo_ASV_combined$sample_names))
match_type_table = matrix(nrow = n_p1_seqs, ncol = n_samples)
count_table_names = c()
unclassified_count_vec = rep(0,n_samples)
g1_unclassified_count_vec = rep(0,n_samples)
for (i in 1:n_p1_seqs) {
  #i=4: "31,30" --> length(p1_seq_split)>1 --> possible_STs==("-","-","8","215") --> Unclassified!
  p1_seq = p1_seqs[i]
  p1_seq_split = strsplit(p1_seq,',')[[1]]
  if (length(p1_seq_split)>1) {
    possible_STs = as.vector(ST_amplicon_table$ST)[which(ST_amplicon_table$epi01_ASV %in% p1_seq_split)]
    if (length(possible_STs)==1) {
      p1_seq = p1_seq_split[1]
    } else {
      p1_seq = "Unclassified"
    }
  }
  if (p1_seq != "Unclassified") {
    p1_seq_ST_table = ST_amplicon_table[which(ST_amplicon_table$epi01_ASV==p1_seq),]
    unique_p1_ASVs = unique(as.vector(p1_seq_ST_table$ST))
    p2_ASVs_matching_p1 = p1_seq_ST_table$epi02_ASV
    count_vec = rep(0,n_samples)
    for (j in 1:n_samples) {
      p1_percent = epidome_object_norm$p1_table[i,j]
      p1_count = eo_ASV_combined$p1_table[i,j]
      if (p1_percent > 0.01) {
        p2_seqs_present_within_difference_threshold_idx = which(epidome_object_norm$p2_table[,j]>(p1_percent-10) & epidome_object_norm$p2_table[,j]>1)
        p2_seqs_present_ASVs = p2_seqs[p2_seqs_present_within_difference_threshold_idx]
        p2_seqs_present_ASVs_matching_p1 = p2_seqs_present_ASVs[which(p2_seqs_present_ASVs %in% p2_ASVs_matching_p1)]
        p2_percent = epidome_object_norm$p2_table[which(p2_seqs %in% p2_seqs_present_ASVs_matching_p1),j]
        p2_count = eo_ASV_combined$p2_table[which(p2_seqs %in% p2_seqs_present_ASVs_matching_p1),j]
        if (length(p2_seqs_present_ASVs_matching_p1) == 0) {
          ST_order = order(p1_seq_ST_table$freq,decreasing = T)
          classification_group = as.vector(p1_seq_ST_table$ST)[ST_order[1]]
          if (classification_group %in% count_table_names) {
            classification_idx = which(count_table_names == classification_group)
            count_table[classification_idx,j] = count_table[classification_idx,j]+p1_count
          } else {
            count_vec = rep(0,n_samples)
            count_vec[j] = p1_count
            count_table = rbind(count_table,count_vec)
            count_table_names = c(count_table_names,classification_group)
          }
          #match_type_table[i,j] = "epi01 match without epi02 match"
          # Replace the match type with actual epi01 and epi02 values
          match_type_table[i, j] = paste0("epi01:", p1_seq, " epi02:NA", " ", classification_group)
        }
        else if (length(p2_seqs_present_ASVs_matching_p1) == 1) {
          p1_p2_seq_ST_table = p1_seq_ST_table[which(p1_seq_ST_table$epi02_ASV==p2_seqs_present_ASVs_matching_p1[1]),]
          ST_order = order(p1_p2_seq_ST_table$freq,decreasing = T)
          classification_group = as.vector(p1_p2_seq_ST_table$ST)[ST_order[1]]
          if (classification_group %in% count_table_names) {
            classification_idx = which(count_table_names == classification_group)
            count_table[classification_idx,j] = count_table[classification_idx,j]+p1_count
          } else {
            count_vec = rep(0,n_samples)
            count_vec[j] = p1_count
            count_table = rbind(count_table,count_vec)
            count_table_names = c(count_table_names,classification_group)
          }
          #match_type_table[i,j] = "Unique epi01 epi02 combination"
          # Replace the match type with actual epi01 and epi02 values
          match_type_table[i, j] = paste0("epi01:", p1_seq, " epi02:", p2_seqs_present_ASVs_matching_p1[1], " ", classification_group)
        } else {
          if (strict_classifier) {
            unclassified_count_vec[j] = unclassified_count_vec[j]+p1_count
          } else {
            p1_p2_seq_ST_table = p1_seq_ST_table[which(p1_seq_ST_table$epi02_ASV %in% p2_seqs_present_ASVs_matching_p1),]
            ST_order = order(p1_p2_seq_ST_table$freq,decreasing = T)
            classification_group = as.vector(p1_p2_seq_ST_table$ST)[ST_order[1]]
            if (classification_group %in% count_table_names) {
              classification_idx = which(count_table_names == classification_group)
              count_table[classification_idx,j] = count_table[classification_idx,j]+p1_count
            } else {
              count_vec = rep(0,n_samples)
              count_vec[j] = p1_count
              count_table = rbind(count_table,count_vec)
              count_table_names = c(count_table_names,classification_group)
            }
          }

          ##unclassified_count_vec[j] = unclassified_count_vec[j]+p1_count
          #match_type_table[i,j] = c("Non unique epi01 epi02 combination")
          # Concatenate all found epi02 values and replace the match type
          epi02_values = paste(p2_seqs_present_ASVs_matching_p1, collapse = ",")
          match_type_table[i, j] = paste0("epi01:", p1_seq, " epi02:", epi02_values, " ", classification_group)
        }

      } else {
        match_type_table[i,j] = "Low counts"
        unclassified_count_vec[j] = unclassified_count_vec[j]+p1_count
      }

    }
  } else {
    count_vec = as.numeric(as.vector(eo_ASV_combined$p1_table[i,]))
    unclassified_count_vec = unclassified_count_vec+count_vec
    g1_unclassified_count_vec = g1_unclassified_count_vec+count_vec
    match_type_vec = rep("Unclassified epi01 match",n_samples)
    match_type_table[i,] = match_type_vec
  }

}
count_table = rbind(count_table,unclassified_count_vec)
rownames(count_table) = c(count_table_names,"Unclassified")
count_table_df <- as.data.frame(count_table)
write.csv(count_table_df, file="count_table_df.csv")
#TODO: delete the last record "Unclassified" due to the few reads!
#count_table_df = count_table_df[-32,]

# --> FOR STEP4: prepare the file match_type_table.csv for treeHeatmap drawing.
write.csv(match_type_table, file = "match_type_table.csv", row.names = TRUE)

count_table_df_ordered = count_table_df[order(rowSums(count_table_df),decreasing = T),]
row_sums <- apply(count_table_df_ordered, MARGIN=1, FUN=sum)
#TODO: the ST in which the counts <= 500 not shown!
#NOTE: ST225 is "epi01:3 epi02:3 225"!
# Since "epi01:3 epi02:NA 225","epi01:3 epi02:NA 225" are 225 disappeared in the match_type_table.csv_____
#Note that epi01 is g216, epi02 is yycH.
#        215            -            8           59          130          297 
#      611637       401179       341187       337409       316389       316273 
#         331            2           73          384          200            5 
#      178616       125734        99899        97885        93613        92385 
#          14          218           23          278          487          290 
#       71417        62467        60192        57213        53630        51454 
#          87           89          100          329          153           83 
#       49000        23390        21880        21635        19453        10539 
#          19          170          225*         570          184           88 
#        9465         1839         1487          863          437          318 
#         114 Unclassified 
#         108           11 
#> row_sums
#         215            -           59 Unclassified          297          130 
#      611915       442192       359708       357868       335880       315734 
#           5           73          331          384          200           14 
#      145578       134301       127734       126465        94917        85440 
#           2          218           23          278          290           87 
#       85438        62492        60192        57213        51363        49000 
#          89          100          329          210           83           19 
#       23390        21880        21635        12663        10642         9465 
#         170          225          570           88          114           10 
#        1839         1487          863          318          108           91 
# > unique(sort(eo_ASV_combined$p1_seqs))
#  [1] "seq1"  "seq10" "seq11" "seq14" "seq16" "seq18" "seq19" "seq2"  "seq20"
# [10] "seq21" "seq22" "seq24" "seq26" "seq28" "seq29" "seq3"  "seq31" "seq37"
# [19] "seq40" "seq5" 
# > unique(sort(eo_ASV_combined$p2_seqs))
#  [1] "seq10" "seq11" "seq12" "seq14" "seq15" "seq18" "seq19" "seq2"  "seq21"
# [10] "seq23" "seq25" "seq26" "seq3"  "seq34" "seq35" "seq36" "seq38" "seq42"
# [19] "seq43" "seq49" "seq5"  "seq52" "seq53" "seq55" "seq56" "seq61" "seq62"
# [28] "seq63" "seq65" "seq7"
#row.names(count_table_df) <- c("ST215","ST130","ST278","ST200","ST184","ST5","ST59","ST83","ST487","ST114","ST8","ST297","ST153","ST384","ST2","ST14","ST89","ST570","-","ST290","ST331","ST73","ST88","ST100","ST87","ST23","ST218","ST329","ST19","ST225","ST170")
#row.names(count_table_df) <- c("215","130","278","200","184","5","59","83","487","114","8","297","153","384","2","14","89","570","-","290","331","73","88","100","87","23","218","329","19","225","170")
#row.names(count_table_df) <- c("ST215","ST130","ST278","ST200","ST5","ST59","ST83","ST114","ST297","ST384","ST14","ST89","ST210","-","ST328","ST331","ST73","ST2","ST88","ST100","ST10","ST290","ST87","ST23","ST218","ST329","ST19","ST225","ST170","Unclassified")

col_order <- c("P7.Nose","P8.Nose","P9.Nose","P10.Nose","P11.Nose","P12.Nose","P13.Nose","P14.Nose","P15.Nose",     "P16.Foot","P17.Foot","P18.Foot","P19.Foot","P20.Foot","P16.Nose","P17.Nose","P18.Nose","P19.Nose","P20.Nose",    "AH.LH","AH.NLH","AH.Nose", "AL.LH","AL.NLH","AL.Nose", "CB.LH","CB.NLH","CB.Nose", "HR.LH","HR.NLH","HR.Nose", "KK.LH","KK.NLH","KK.Nose", "MC.LH","MC.NLH","MC.Nose", "MR.LH","MR.NLH","MR.Nose",  "PT2.LH","PT2.NLH","PT2.Nose", "RP.LH","RP.NLH","RP.Nose", "SA.LH","SA.NLH","SA.Nose", "XN.LH","XN.NLH","XN.Nose") 
count_table_df_reordered <- count_table_df_ordered[,col_order]
p = make_barplot_epidome(count_table_df_reordered,reorder=FALSE,normalize=TRUE)
#p = make_barplot_epidome(count_table_df_reordered,reorder=TRUE,normalize=TRUE)
png("Barplot_All.png", width=1600, height=900)
p
dev.off()

#Change rowname from '-' to 'Novel'
rownames(count_table_df_reordered)[rownames(count_table_df_reordered) == "-"] <- "Novel"
write.csv(file="count_table_df_reordered.txt", count_table_df_reordered)

  4. prepare data for plotTreeHeatmap: table+tree.

      #replace "," to \n in match_type_table.csv
  sed 's/\",\"/\n/g' match_type_table.csv > match_type_table.csv_
  grep "epi02" match_type_table.csv_ > match_type_table.csv__
  grep -v "NA" match_type_table.csv__ > match_type_table.csv___
  #(Deprecated) awk -F' ' '{ split($2, a, ":"); split(a[2], b, ","); print $1 " epi02:" b[1] " " $3 }' match_type_table.csv___ > match_type_table.csv____
  #(Deprecated) sed -i -e 's/\"//g' match_type_table.csv____
  #(Deprecated) sort match_type_table.csv____ -u > match_type_table.csv_____

  #Manually select and merge items from match_type_table.csv___ and save it to match_type_table.csv____ (see content as below)!
      epi01:1 epi02:36 2
      epi01:3 epi02:3 225
      epi01:5 epi02:34 130
      epi01:5 epi02:3 278
      epi01:5 epi02:43 200
      epi01:5 epi02:36 184
      epi01:11 epi02:19 19
      epi01:14 epi02:34 297
      epi01:14 epi02:36 153
      epi01:16 epi02:56 329
      epi01:16 epi02:55 329
      epi01:18 epi02:14 218
      epi01:19 epi02:43 170
      epi01:2 epi02:15 -
      epi01:20 epi02:36 2
      epi01:20 epi02:43 14
      epi01:20 epi02:42 89
      epi01:20 epi02:34 384
      epi01:20 epi02:11 570
      epi01:21 epi02:34 100
      epi01:21 epi02:63 290
      epi01:24 epi02:38 23
      epi01:24 epi02:2 87
      epi01:26 epi02:34 -
      epi01:26 epi02:53 290
      epi01:26 epi02:63 331
      epi01:28 epi02:43 215
      epi01:29 epi02:43 215
      epi01:29 epi02:55 -
      epi01:31 epi02:43 8
      epi01:37 epi02:36 2
      epi01:37 epi02:26 88
      epi01:37 epi02:43 73
      epi01:40 epi02:2 5
      epi01:40 epi02:34 5
      epi01:40 epi02:34 59
      epi01:40 epi02:26 114
      epi01:40 epi02:43 83
      epi01:40 epi02:36 487
      #---- v2 ('-' three times, ST215, ST290, ST329, ST59 two times) --> delete 'epi01:28 epi02:43 215' ----
      epi01:1 epi02:34 -
      epi01:3 epi02:3 225
      epi01:5 epi02:3 278
      epi01:5 epi02:34 130
      epi01:5 epi02:43 200
      epi01:11 epi02:19 19
      epi01:14 epi02:34 297
      epi01:16 epi02:55 329
      epi01:16 epi02:56 329
      epi01:18 epi02:14 218
      epi01:19 epi02:43 170
      epi01:20 epi02:42 89
      epi01:20 epi02:43 14
      epi01:20 epi02:34 384
      epi01:20 epi02:11 570
      epi01:20 epi02:3 210
      epi01:21 epi02:38 10
      epi01:21 epi02:34 100
      epi01:21 epi02:63 290
      epi01:24 epi02:38 23
      epi01:24 epi02:2 87
      epi01:26 epi02:34 -
      epi01:26 epi02:53 290
      epi01:26 epi02:63 331
      epi01:29 epi02:43 215
      epi01:29 epi02:55 -
      epi01:37 epi02:62 2
      epi01:37 epi02:26 88
      epi01:37 epi02:43 73
      epi01:37 epi02:34 59
      epi01:40 epi02:2 5
      epi01:40 epi02:26 114
      epi01:40 epi02:43 83
      epi01:40 epi02:34 59

python3 concat_epi01_epi02.py
    # Define the function to extract fasta sequences
    def fasta_to_dict(filename):
        sequences = {}
        header = None
        sequence = []
        with open(filename, 'r') as file:
            for line in file:
                line = line.strip()
                if line.startswith('>'):
                    if header:
                        sequences[header] = ''.join(sequence)
                    header = line[1:]
                    sequence = []
                else:
                    sequence.append(line)
            if header:
                sequences[header] = ''.join(sequence)
        return sequences
    # Read combinations from match_type_table.csv____
    combinations = []
    with open("match_type_table.csv____", 'r') as csv_file:
        for line in csv_file:
            parts = line.strip().split()
            epi01 = parts[0].split(":")[1]
            epi02 = parts[1].split(":")[1]
            st = parts[2]
            combinations.append((epi01, epi02, st))

    # Read both fasta files into dictionaries
    g216_sequences = fasta_to_dict("../epidome/DB/epi01_ref_aln.fasta")
    yycH_sequences = fasta_to_dict("../epidome/DB/epi02_ref_aln.fasta")

    output_filename = "concatenated_sequences.fasta"
    with open(output_filename, 'w') as output_file:
        for epi01, epi02, st in combinations:
            if epi01 in g216_sequences and epi02 in yycH_sequences:
                concatenated_seq = g216_sequences[epi01] + "NNN" + yycH_sequences[epi02]
                header = f">g216-{epi01}_yycH-{epi02}_ST{st}"
                output_file.write(f"{header}\n{concatenated_seq}\n")

    print(f"Concatenated sequences saved to {output_filename}")

sed -i -e 's/ST-/-/g' concatenated_sequences.fasta
muscle -in concatenated_sequences.fasta -out aligned.fasta
muscle -clw -in concatenated_sequences.fasta -out aligned.clw
FastTree -nt < aligned.fasta > concatenated_sequences.tree
# Run plotTree.R to draw ggtree_and_gheatmap.png.

library(ggtree)
library(ggplot2)

setwd("/home/jhuang/DATA/Data_Luise_Epidome_batch3/run_2023_98_2/")

# -- edit tree --
#https://icytree.org/
info <- read.csv("typing_34.csv")
info$name <- info$Isolate
tree <- read.tree("concatenated_sequences.tree")
#cols <- c(infection='purple2', commensalism='skyblue2')     

library(dplyr)
heatmapData2 <- info %>% select(Isolate,  g216, ST)
rn <- heatmapData2$Isolate
heatmapData2$Isolate <- NULL
heatmapData2 <- as.data.frame(sapply(heatmapData2, as.character))
rownames(heatmapData2) <- rn

#https://bookdown.org/hneth/ds4psy/D-3-apx-colors-basics.html
#"g216:16","g216:18","g216:19","g216:20"(4),"g216:21"(2),"g216:24"(1),"g216:26"(2)
heatmap.colours <- c("darkgrey",   "palegreen","seagreen","seagreen","seagreen4","seagreen4","yellowgreen",    "orange4","orange4","orange4","orange4","orange4",    "dodgerblue3","blue4","dodgerblue3","dodgerblue3",  "azure4",    "magenta2","maroon2","mediumorchid1","mediumorchid1","mediumorchid1","mediumorchid1","mediumorchid1", "mediumorchid3","mediumorchid3","mediumorchid3","mediumpurple1","mediumpurple4","mediumpurple4",         "red4",                                            "olivedrab3","olivedrab4","seagreen","seagreen4","yellowgreen",    "orange4",    "blue4","dodgerblue3",           "magenta2","maroon2","mediumorchid1", "mediumorchid3","mediumpurple1","mediumpurple4",         "red4")
names(heatmap.colours) <- c("-",    "ST225","ST130","ST278","ST19","ST200", "ST297",                                             "ST5","ST59","ST114","ST83","ST487",   "ST2","ST215","ST73","ST88",   "ST8",                                                               "ST329","ST170", "ST89","ST14","ST570","ST384","ST210",                  "ST10","ST100","ST290","ST23","ST87", "ST331",     "ST218",      "g216:1","g216:3","g216:5","g216:11","g216:14",   "g216:40",   "g216:29","g216:37",        "g216:16","g216:19","g216:20","g216:21","g216:24","g216:26",    "g216:18")
#mydat$Regulation <- factor(mydat$Regulation, levels=c("up","down"))

#circular
#geom_tippoint(aes(color=Type)) + scale_color_manual(values=cols) + 
p <- ggtree(tree, layout='circular', branch.length='none') %<+% info + geom_tiplab2(aes(label=name), offset=1)
#, geom='text', align=TRUE,  linetype=NA, hjust=1.8,check.overlap=TRUE, size=3.3
#difference between geom_tiplab and geom_tiplab2?
#+ theme(axis.text.x = element_text(angle = 30, vjust = 0.5)) + theme(axis.text = element_text(size = 20))  + scale_size(range = c(1, 20))
#font.size=10, 
png("ggtree.png", width=1260, height=1260)
#svg("ggtree.svg", width=1260, height=1260)
p
dev.off()

desired_order <- c("-",    "ST225","ST130","ST278","ST19","ST200", "ST297",                                             "ST5","ST59","ST114","ST83","ST487",   "ST2","ST215","ST73","ST88",   "ST8",                                                               "ST329","ST170", "ST89","ST14","ST570","ST384","ST210",                  "ST10","ST100","ST290","ST23","ST87", "ST331",     "ST218",      "g216:1","g216:3","g216:5","g216:11","g216:14",   "g216:40",   "g216:29","g216:37",        "g216:16","g216:19","g216:20","g216:21","g216:24","g216:26",    "g216:18")

png("ggtree_and_gheatmap.png", width=1290, height=1000)
#svg("ggtree_and_gheatmap.svg", width=17, height=15)
gheatmap(p, heatmapData2, width=0.15,colnames_position="top", colnames_angle=90, colnames_offset_y = 0.1, hjust=0.5, font.size=4, offset = 9) + scale_fill_manual(values=heatmap.colours, breaks=desired_order) +  theme(legend.text = element_text(size = 14)) + theme(legend.title = element_text(size = 14)) + guides(fill=guide_legend(title=""), color = guide_legend(override.aes = list(size = 5)))  
dev.off()

  5. prepre the table for the alpha diversity calculation

    To calculate the alpha diversity values for each row based on the provided metrics (chao1, observed_otus, shannon, PD_whole_tree), you’ll typically use specialized software packages such as QIIME or R packages like vegan. Here, I can provide a high-level approach and calculations for some of the metrics. However, some of the metrics like PD_whole_tree require phylogenetic trees, which can’t be derived from the table you provided.

    import pandas as pd
import numpy as np

# Read the data from CSV
df = pd.read_csv('count_table_df_reordered.csv', index_col=0)

def chao1(s):
s = np.array(s)
n = np.sum(s)
s_obs = np.count_nonzero(s)
s1 = np.sum(s == 1)
s2 = np.sum(s == 2)
return s_obs + ((s1 ** 2) / (2 * s2)) if s2 != 0 else s_obs

def observed_sts(s):
return np.count_nonzero(s)

def shannon2(s):
s = np.array(s)
p = s / np.sum(s)
return -np.sum(p * np.log(p + 1e-10))  # added small value to avoid log(0)

def shannon(s):
s = np.array(s)
total = np.sum(s)
p = s / total  # Convert absolute counts to proportions
return -np.sum(p * np.log(p + 1e-10))  # added small value to avoid log(0)

# Calculate the metrics for each sample (column)
results = {
'sample': [],
'chao1': [],
'observed_sts': [],
'shannon': []
}

for column in df.columns:
results['sample'].append(column)
results['chao1'].append(chao1(df[column]))
results['observed_sts'].append(observed_sts(df[column]))
results['shannon'].append(shannon(df[column]))

# Convert results to DataFrame and save to a new CSV
result_df = pd.DataFrame(results)
result_df.to_csv('alpha_diversity_metrics_samples.csv', index=False)

    For PD_whole_tree, you would require a phylogenetic tree of the OTUs and then compute the sum of branch lengths. This is more complex and requires specialized software and the phylogenetic tree information, which isn’t present in the table you provided.

    The input data being absolute counts versus relative abundance does influence certain alpha diversity metrics. Let’s re-evaluate:

    • Observed OTUs: Whether the data is in absolute counts or relative abundance doesn’t impact the “Observed OTUs” metric. It’s simply counting how many OTUs have a non-zero count.

    • Chao1: Chao1 is intended for absolute counts. The metric is based on the number of singletons and doubletons, which are derived from raw counts and not relative abundances. Hence, the calculation remains accurate.

    • Shannon Diversity Index: The Shannon Diversity Index is based on the relative abundances of OTUs, not their absolute counts. The formula requires the proportion of each species (or OTU) in the sample. Thus, you’d need to convert absolute counts to proportions (or relative abundance) before calculating this metric.

  6. calculate the alpha diversity in Phyloseq.Rmd. The code for the calculation is as follows.

    \pagebreak

# Alpha diversity for ST
Plot Shannon index and observed STs. 
Regroup together samples from the same group.

```{r, echo=TRUE, warning=FALSE}
hmp.div_st <- read.csv("alpha_diversity_metrics_samples.csv", sep=",") 
colnames(hmp.div_st) <- c("sample", "chao1", "observed_sts", "shannon")
row.names(hmp.div_st) <- hmp.div_st$sample
div.df <- merge(hmp.div_st, hmp.meta, by.x="sample", by.y = "sam_name")
div.df2 <- div.df[, c("Group", "shannon", "observed_sts")]
colnames(div.df2) <- c("Group", "Shannon", "ST")
#colnames(div.df2)
options(max.print=999999)

stat.test.Shannon <- compare_means(
Shannon ~ Group, data = div.df2,
method = "t.test"
)
knitr::kable(stat.test.Shannon) %>% kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive"))

stat.test.ST <- compare_means(
ST ~ Group, data = div.df2,
method = "t.test"
)
knitr::kable(stat.test.ST) %>% kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive"))

div_df_melt <- reshape2::melt(div.df2)

p <- ggboxplot(div_df_melt, x = "Group", y = "value",
              facet.by = "variable", 
              scales = "free",
              width = 0.5,
              fill = "gray", legend= "right")
#ggpar(p, xlab = FALSE, ylab = FALSE)
lev <- levels(factor(div_df_melt$Group)) # get the variables
L.pairs <- combn(seq_along(lev), 2, simplify = FALSE, FUN = function(i) lev[i]) #%>% filter(p.signif != "ns")
my_stat_compare_means  <- function (mapping = NULL, data = NULL, method = NULL, paired = FALSE, 
    method.args = list(), ref.group = NULL, comparisons = NULL, 
    hide.ns = FALSE, label.sep = ", ", label = NULL, label.x.npc = "left", 
    label.y.npc = "top", label.x = NULL, label.y = NULL, tip.length = 0.03, 
    symnum.args = list(), geom = "text", position = "identity", 
    na.rm = FALSE, show.legend = NA, inherit.aes = TRUE, ...) 
{
    if (!is.null(comparisons)) {
        method.info <- ggpubr:::.method_info(method)
        method <- method.info$method
        method.args <- ggpubr:::.add_item(method.args, paired = paired)
        if (method == "wilcox.test") 
            method.args$exact <- FALSE
        pms <- list(...)
        size <- ifelse(is.null(pms$size), 0.3, pms$size)
        color <- ifelse(is.null(pms$color), "black", pms$color)
        map_signif_level <- FALSE
        if (is.null(label)) 
            label <- "p.format"
        if (ggpubr:::.is_p.signif_in_mapping(mapping) | (label %in% "p.signif")) {
            if (ggpubr:::.is_empty(symnum.args)) {
                map_signif_level <- c(`****` = 1e-04, `***` = 0.001, 
                  `**` = 0.01, `*` = 0.05, ns = 1)
            } else {
              map_signif_level <- symnum.args
            } 
            if (hide.ns) 
                names(map_signif_level)[5] <- " "
        }
        step_increase <- ifelse(is.null(label.y), 0.12, 0)
        ggsignif::geom_signif(comparisons = comparisons, y_position = label.y, 
            test = method, test.args = method.args, step_increase = step_increase, 
            size = size, color = color, map_signif_level = map_signif_level, 
            tip_length = tip.length, data = data)
    } else {
        mapping <- ggpubr:::.update_mapping(mapping, label)
        layer(stat = StatCompareMeans, data = data, mapping = mapping, 
            geom = geom, position = position, show.legend = show.legend, 
            inherit.aes = inherit.aes, params = list(label.x.npc = label.x.npc, 
                label.y.npc = label.y.npc, label.x = label.x, 
                label.y = label.y, label.sep = label.sep, method = method, 
                method.args = method.args, paired = paired, ref.group = ref.group, 
                symnum.args = symnum.args, hide.ns = hide.ns, 
                na.rm = na.rm, ...))
    }
}

p3 <- p + 
  stat_compare_means(
    method="t.test",
    comparisons = list(c("P16-P20.Foot", "P16-P20.Nose"), c("AH-XN.LH", "AH-XN.Nose"), c("AH-XN.NLH", "AH-XN.Nose")), 
    label = "p.signif",
    symnum.args <- list(cutpoints = c(0, 0.0001, 0.001, 0.01, 0.05, 1), symbols = c("****", "***", "**", "*", "ns")),
  ) +
  theme(axis.text.x = element_text(angle = 30, hjust = 1))  # add this line to rotate x-axis text
print(p3)

ggsave("./figures/alpha_diversity_ST_Group2.png", device="png", height = 10, width = 12)
ggsave("./figures/alpha_diversity_ST_Group2.svg", device="svg", height = 10, width = 12)
```

Leave a Reply

Your email address will not be published. Required fields are marked *