<jupyter_start><jupyter_text>Table of Contents1&nbsp;&nbsp;Goal2&nbsp;&nbsp;Var3&nbsp;&nbsp;Init4&nbsp;&nbsp;Load5&nbsp;&nbsp;Compare pass/fail5.1&nbsp;&nbsp;check/format5.2&nbsp;&nbsp;Summary5.2.1&nbsp;&nbsp;Success rate5.2.1.1&nbsp;&nbsp;Writing table of samples5.2.2&nbsp;&nbsp;No. of samples5.2.3&nbsp;&nbsp;Passed per species5.2.4&nbsp;&nbsp;Mapping onto host tree5.2.4.1&nbsp;&nbsp;Host metadata6&nbsp;&nbsp;sessionInfo  Goal* Summary of NGS pass/fail samples  * how many succeeded of the total?    Var<jupyter_code>work_dir = '/ebio/abt3_projects/Georg_animal_feces/data/16S_arch/NGS/'

# all samples
metadata_file = '/ebio/abt3_projects/Georg_animal_feces/data/mapping/unified_metadata_complete_200429.tsv'

# host tree for all species
host_tree_file = '/ebio/abt3_projects/Georg_animal_feces/data/animal/phylogeny/TT/unified_metadata_complete_20190416_rn.nwk'

## all samples & taxa
physeq_all_file = file.path('/ebio/abt3_projects/Georg_animal_feces/data/16S_arch/MiSeq-Runs-116-122-126-189-190/LLA/merged/phyloseq/physeq_all.RDS')

## 16S-arch phyloseq object (one sample per individual)
physeq_arch_IndD_file = '/ebio/abt3_projects/Georg_animal_feces/data/16S_arch/MiSeq-Runs-116-122-126-189-190/LLA/merged/phyloseq/physeq_arch_IndD.RDS'

# params
threads = 8<jupyter_output><empty_output><jupyter_text>Init<jupyter_code>library(dplyr)
library(tidyr)
library(ggplot2)
library(ape)
library(phyloseq)
library(LeyLabRMisc)
df.dims()<jupyter_output><empty_output><jupyter_text>Load<jupyter_code># all sample metadata
metadata = read.delim(metadata_file, sep='\t')
metadata 
# host tree
host_tree = read.tree(host_tree_file)
host_tree
# all samples
physeq_all = readRDS(physeq_all_file)
## filtering
x = physeq_all %>% sample_names
x = x[grepl('Neg_control', x)]
physeq_all_neg = subset_samples(physeq_all, SampleID %in% x) %>%
    filter_taxa(function(x) max(x) > 0, TRUE)
## summary
physeq_all_neg
# all passed samples
physeq_arch_IndD = readRDS(physeq_arch_IndD_file)
physeq_arch_IndD
# formatting
otu = physeq_arch_IndD %>%
    phyloseq2df(otu_table, long=TRUE) %>%
    group_by(Sample) %>%
    mutate(Total_count = sum(Count)) %>%
    ungroup() %>%
    filter(Count > 0)

tax = physeq_arch_IndD %>%
    phyloseq2df(tax_table) %>%
    filter(Class == 'Bathyarchaeia')

otu = otu %>%
    inner_join(tax, c('OTU')) 

otu %>% .$Sample %>% table
# summary
physeq_arch_IndD %>%
    phyloseq2df(tax_table) %>% 
    distinct(Phylum, Class) %>%
    arrange(Phylum, Class)
# MiSeq runs 116-122-126
samps = c('1_Eurasian_Beaver', '2_Raccoon_Dog', '5_European_Hare', '21_Yellow_necked_Field_Mouse',
  '43_European_Badger', '66_Mouflon', '73_One_humped_Camel', '77_Onager', '78_Indian_Gazelle',
  '81_Wolf', '87_Garden_Dormouse', '88_European_Rabbit', '101_Horse', '102_Kulan', '111_West_European_Hedgehog',
  '122_Wild_Boar', '137_Gaur', '140_Sambar', '144_Cattle', '152_Tree_Shrew', '198_Human',
  '206_Red_Sheep', '237_Red_shanked_Douc_Langur', '272_Black_Grouse', '275_Barnacle_Goose',
  '279_Western_Grey_Kangaroo', '285_Short_beaked_Echidna', '287_Common_Brushtail', '289_Eastern_Grey_Kangaroo',
  '311_Domestic_Cat', '320_Red_Deer', '376_Kulan', '395_Przewalski_horse', '407_African_Bush_Elephant',
  '417_Giraffe', '419_African_Buffalo', '422_Common_Hippopotamus', '424_Quagga', '95_Meadow_Viper',
  '123_Alpine_Chamois', '128_Alpine_Marmot', '129_Alpine_Ibex', '131_Rock_Ptarmigan', '155_Long_tailed_Field_Mouse',
  '169_Great_Cormorant', '268_Rook', '351_Grass_Snake', '389_Indian_Rhinoceros', '7_European_Otter',
  '9_Indian_Gazelle', '11_Onager', '13_Wolf', '15_Tawny_Owl', '16_Fat_Dormouse', '22_European_Rabbit',
  '23_Garden_Dormouse', '31_Chicken', '32_Greylag_Goose', '33_Wild_Turkey', '34_Mallard_Duck',
  '40_Bactrian_Camel', '41_Pika', '42_Goitered_Gazelle', '60_European_Chub', '85_Domestic_Dog', 
  '93_Red_Sheep', '94_Mangalica', '96_European_Hare', '116_Common_Kestrel', '121_Raccoon_Dog', 
  '135_Horse', '138_Gaur', '139_Sambar', '141_Red_cheeked_Gibbon', '150_Cattle', '154_Tree_Shrew', 
  '157_Grey_Heron', '158_White_tailed_Eagle', '164_Greylag_Goose', '170_Tawny_Owl', '172_European_Greenfinch', 
  '174_Great_Tit', '195_Domestic_Cat', '196_Domestic_Cat', '197_Human', '199_Human', '200_Human', 
  '203_Red_Deer', '205_Red_Sheep', '207_Red_Sheep', '208_Red_Sheep', '209_Red_Sheep', '210_Red_Sheep', 
  '211_Red_Sheep', '212_Goose', '213_Goose', '214_Goose', '215_Goose', '218_Goose', '219_Goose', 
  '221_Alpine_Ibex', '223_Human', '224_Domestic_Cat', '225_Domestic_Cat', '226_Domestic_Cat',
  '228_Domestic_Cat', '230_Carrion_Crow', '234_Beech_Marten', '238_Hanuman_Langur', 
  '240_Southern_White_cheeked_Gibbon', '252_Eurasian_Lynx', '259_Mute_Swan', '263_White_Stork', 
  '265_Western_Marsh_Harrier', '270_Common_Pheasant', '271_Common_Pheasant', '276_Barnacle_Goose', 
  '280_Western_Grey_Kangaroo', '283_Koala', '284_Koala', '286_Short_beaked_Echidna', '290_Eastern_Grey_Kangaroo', 
  '301_European_Hare', '302_European_Hare', '303_European_Hare', '304_European_Hare', '305_European_Hare', 
  '306_European_Hare', '307_European_Hare', '308_European_Hare', '309_Domestic_Cat', '310_Domestic_Cat', 
  '312_Domestic_Cat', '313_Red_Deer', '314_Red_Deer', '315_Red_Deer', '316_Red_Deer', '317_Red_Deer', 
  '318_Red_Deer', '319_Red_Deer', '322_Red_Deer', '324_Red_Deer', '330_Ide', '332_Harbor_Porpoise', 
  '333_Gray_Seal', '336_Sei_Whale', '338_White_beaked_Dolphin', '339_White_beaked_Dolphin', 
  '340_Western_Lowland_Gorilla', '341_Western_Lowland_Gorilla', '344_Indian_Rhinoceros', '347_Brown_Bear', 
  '348_European_Rabbit', '352_Kinkajou', '353_Raccoon', '355_Hoffmanns_Two_toed_Sloth', '356_Northern_Tamandua', 
  '360_White_tailed_Deer', '362_Gray_Four_eyed_Opossum', '363_Desmarests_Spiny_Pocket_Mouse', 
  '368_European_Ground_Squirrel', '369_Koala', '370_African_Bush_Elephant', '371_Linnaeus_Two_toed_Sloth', 
  '372_Emu', '378_Fat_Sand_Rat', '379_Koala', '390_European_Ground_Squirrel', '391_European_Ground_Squirrel', 
  '393_European_Rabbit', '394_Przewalski_horse', '396_Morelets_crocodile', '397_Penguin', '404_Brown_Greater_Galago', 
  '408_African_Bush_Elephant', '409_Impala', '410_Impala', '411_Leopard', '412_Lion', '414_Vervet_Monkey', 
  '415_White_Rhinoceros', '416_White_Rhinoceros', '418_Giraffe', '420_Cattle', '421_Common_Hippopotamus', 
  '423_Quagga', '426_Ural_Owl', '427_Great_Spotted_Woodpecker', '428_African_Bush_Elephant', '429_Giraffe', 
  '430_Giraffe', '432_Vervet_Monkey', '433_Raccoon', '434_Blue_Wildebeest', '435_Bactrian_Camel', 
  'F14_Common_Bream', 'F35_Red_Deer', 'F36_Red_Deer', 'F44_Fallow_Deer', 'F45_Red_Deer', 'F46_Red_Deer', 
  'F47_Red_Deer', 'F48_Red_Deer', 'Pos_control_7', 'Neg_control_7', 'F53_Mouflon', 'F66_Wild_Boar', 
  'F68_Red_Deer', 'F69_Red_Deer', 'F70_Red_Deer', 'F80_Red_Deer', 'F90_Domestic_Dog', '108_Striped_Desert_Hamster', 
  '109_Red_Fox', '117_Long_eared_Owl', '119_Ural_Owl', '192_Common_Frog', '194_Agile_Frog', '125_Common_Carp', 
  '126_White_Stork', '130_Alpine_Ibex', '146_Cattle', '153_Bank_Vole', '179_Roach', '233_Aesculapian_Snake', 
  '236_Pygmy_Slow_Loris', '241_Italian_wall_lizard', '242_Dalmatian_Tortoise', '243_Noctule', 
  '247_Common_Hamster', '260_Blackbird', '264_Eurasian_Magpie', '266_Western_Marsh_Harrier', 
  '288_Common_Brushtail', '292_Central_Bearded_Dragon', '297_Horsfields_Bronze_Cuckoo', 
  '326_Silver_Bream', '327_Silver_Bream', '328_Prussian_Carp', '349_Smooth_Newt', '350_Smooth_Newt', 
  '380_Wild_Boar', '381_Wild_Boar', '382_Arctic_Wolf', '383_Arctic_Wolf', 'Pos_control_9', 
  '384_Eurasian_Lynx', '386_Scalloped_Hammerhead', '400_House_Mouse', 'Neg_control_9', '43_European_Badger', 
  '66_Mouflon', '73_One_humped_Camel', '77_Onager', '88_European_Rabbit', '101_Horse', '102_Kulan', 
  '111_West_European_Hedgehog', '122_Wild_Boar', '137_Gaur', '140_Sambar', '144_Cattle', '206_Red_Sheep', 
  '237_Red_shanked_Douc_Langur', '279_Western_Grey_Kangaroo', '285_Short_beaked_Echidna', '287_Common_Brushtail', 
  '289_Eastern_Grey_Kangaroo', '320_Red_Deer', '376_Kulan', '395_Przewalski_horse', '407_African_Bush_Elephant', 
  '417_Giraffe', '419_African_Buffalo', '422_Common_Hippopotamus', '424_Quagga', '95_Meadow_Viper', 
  '123_Alpine_Chamois', '128_Alpine_Marmot', '129_Alpine_Ibex', '131_Rock_Ptarmigan', '268_Rook', 
  '351_Grass_Snake', '389_Indian_Rhinoceros', '7_European_Otter', '11_Onager', '13_Wolf', '23_Garden_Dormouse', 
  '31_Chicken', '32_Greylag_Goose', '34_Mallard_Duck', '40_Bactrian_Camel', '42_Goitered_Gazelle', 
  '60_European_Chub', '94_Mangalica', '96_European_Hare', '116_Common_Kestrel', '121_Raccoon_Dog', 
  '135_Horse', '138_Gaur', '139_Sambar', '150_Cattle', '154_Tree_Shrew', '172_European_Greenfinch', 
  '174_Great_Tit', '203_Red_Deer', '205_Red_Sheep', '207_Red_Sheep', '208_Red_Sheep', '209_Red_Sheep', 
  '210_Red_Sheep', '211_Red_Sheep', '212_Goose', '213_Goose', '214_Goose', '215_Goose', '218_Goose', 
  '219_Goose', '221_Alpine_Ibex', '238_Hanuman_Langur', '240_Southern_White_cheeked_Gibbon', 
  '270_Common_Pheasant', '271_Common_Pheasant', '280_Western_Grey_Kangaroo', '283_Koala', 
  '286_Short_beaked_Echidna', '290_Eastern_Grey_Kangaroo', '313_Red_Deer', '314_Red_Deer', 
  '315_Red_Deer', '316_Red_Deer', '317_Red_Deer', '318_Red_Deer', '319_Red_Deer', '322_Red_Deer', 
  '324_Red_Deer', '330_Ide', '333_Gray_Seal', '336_Sei_Whale', '339_White_beaked_Dolphin', 
  '340_Western_Lowland_Gorilla', '341_Western_Lowland_Gorilla', '344_Indian_Rhinoceros', 
  '352_Kinkajou', '356_Northern_Tamandua', '360_White_tailed_Deer', '368_European_Ground_Squirrel', 
  '369_Koala', '370_African_Bush_Elephant', '371_Linnaeus_Two_toed_Sloth', '372_Emu', '378_Fat_Sand_Rat', 
  '390_European_Ground_Squirrel', '396_Morelets_crocodile', '408_African_Bush_Elephant', '409_Impala', 
  '410_Impala', '411_Leopard', '415_White_Rhinoceros', '416_White_Rhinoceros', '418_Giraffe', '420_Cattle', 
  '421_Common_Hippopotamus', '423_Quagga', '428_African_Bush_Elephant', '429_Giraffe', '430_Giraffe', 
  '432_Vervet_Monkey', '434_Blue_Wildebeest', '435_Bactrian_Camel', 'F14_Common_Bream', 'F35_Red_Deer', 
  'F36_Red_Deer', 'F44_Fallow_Deer', 'F45_Red_Deer', 'F46_Red_Deer', 'F47_Red_Deer', 'F48_Red_Deer', 
  'F53_Mouflon', 'F66_Wild_Boar', 'F68_Red_Deer', 'F69_Red_Deer', 'F70_Red_Deer', 'F80_Red_Deer', 
  'F90_Domestic_Dog', '117_Long_eared_Owl', '130_Alpine_Ibex', '146_Cattle', '179_Roach', 
  '233_Aesculapian_Snake', '236_Pygmy_Slow_Loris', '242_Dalmatian_Tortoise', '247_Common_Hamster', 
  '266_Western_Marsh_Harrier', '288_Common_Brushtail', '327_Silver_Bream', '328_Prussian_Carp', 
  '349_Smooth_Newt', '350_Smooth_Newt', '380_Wild_Boar', '381_Wild_Boar', '382_Arctic_Wolf')

samps = unique(samps)
samps = samps[!grepl('_control_', samps)]
samps %>% length
# function to read in PCR results
read_file = function(F, D){
  x = read.delim(file.path(D, F), sep='\t')
  x$labware_name = gsub('.+(plate[0-9]).+', '\\1', F)
  x$labware_type = '96 Well Eppendorf TwinTec PCR'
  return(x)
}
# PCR results
rfu_files = c(
  'sd_2019-11-05_11-00-40_Connect1_1.PCR_Georg_16S_Archaea_plate1_2.part_End_Point_Results.txt',
  'sd_2019-11-05_11-03-35_Connect2_1.PCR_Georg_16S_Archaea_plate2_2.part_End_Point_Results.txt',
  'sd_2019-11-05_11-06-07_Connect3_1.PCR_Georg_16S_Archaea_plate3_2.part_End_Point_Results.txt'
)
D = file.path(work_dir, '191000_SilkeD-prep', 'PCR-step1_plate2')
rfu_sd = rfu_files %>%
  as.list %>%
  lapply(read_file, D=D) %>%
  do.call(rbind, .) %>%
  filter(!grepl('control', Sample)) 
rfu_sd
# plate maps
F = file.path(work_dir, '191000_SilkeD-prep', 'PCR_step1_plate1', 'TECAN_NGS_amplicon_PCR1_map.txt')
plate_map = read.delim(F, sep='\t') %>%
    filter(!grepl('control', SampleID)) %>%
    mutate(SampleID = SampleID %>% as.character)
plate_map
# summary
all_samps = c(as.character(rfu_sd$Sample), samps, plate_map$SampleID)
all_samps = all_samps[!grepl('control', all_samps)]
all_samps = gsub('^([^X])', 'X\\1', all_samps)
all_samps = gsub('^(^XF)', 'F', all_samps)
all_samps = all_samps %>% unique %>% sort
all_samps %>% length %>% print
all_samps %>% print<jupyter_output>[1] 313
  [1] "F14_Common_Bream"                   "F141_Wood_Sandpiper"               
  [3] "F157a_European_Toad"                "F157b_European_Toad"               
  [5] "F241_Great_Cormorant"               "F300_Asp"                          
  [7] "F302_Northern_Pike"                 "F35_Red_Deer"                      
  [9] "F36_Red_Deer"                       "F42_Fallow_Deer"                   
 [11] "F44_Fallow_Deer"                    "F45_Red_Deer"                      
 [13] "F46_Red_Deer"                       "F47_Red_Deer"                      
 [15] "F48_Red_Deer"                       "F53_Mouflon"                       
 [17] "F66_Wild_Boar"                      "F68_Red_Deer"                      
 [19] "F69_Red_Deer"                       "F70_Red_Deer"                      
 [21] "F80_Red_Deer"                       "F90_Domestic_Dog"                  
 [23] "X1_Eurasian_Beaver"                 "X10_House_Mouse"                   
 [25] "X101_Horse"              [...]<jupyter_text>Compare pass/fail<jupyter_code>physeq_arch_IndD %>%
    phyloseq2df(otu_table, long = TRUE) %>%
    group_by(Sample) %>%
    summarize(Count = sum(Count), .groups='drop') %>%
    .$Count %>% summary_x<jupyter_output><empty_output><jupyter_text>check/format<jupyter_code># checking overlap
overlap(physeq_arch_IndD %>% sample_names, all_samps)
overlap(metadata$SampleID, all_samps)
overlap(metadata$SampleID, all_samps, to_return='diff_y')
metadata = metadata %>%
    filter(SampleID %in% all_samps) %>%
    mutate(NGS_pass = SampleID %in% (physeq_arch_IndD %>% sample_names))
metadata
# summary
metadata$NGS_pass %>% summary
# filter by class
metadata = metadata %>%
    filter(class != 'Chondrichthyes')
metadata$NGS_pass %>% summary<jupyter_output><empty_output><jupyter_text>Summary  Success rate<jupyter_code>df.dims(20)
metadata %>% 
    group_by(class) %>%
    summarize(n_samples = n(), .groups='drop')
df.dims()
# overall success rate
sr_class = metadata %>%
    group_by(class) %>%
    mutate(n_samples = n()) %>%
    group_by(class, NGS_pass) %>%
    summarize(rate = n() / first(n_samples) * 100, 
              n_samples = first(n_samples), 
              .groups='drop') %>%
    filter(NGS_pass == TRUE) %>%
    mutate(category = 'Class') %>%
    rename('feature' = class)
sr_class
# overall success rate
sr_st = metadata %>%
    group_by(sample_type) %>%
    mutate(n_samples = n()) %>%
    group_by(sample_type, NGS_pass) %>%
    summarize(rate = n() / first(n_samples) * 100, 
              n_samples = first(n_samples), 
              .groups='drop') %>%
    filter(NGS_pass == TRUE) %>%
    mutate(category = 'Sample type') %>%
    rename('feature' = sample_type)
sr_st
# overall success rate
sr_diet = metadata %>%
    group_by(diet) %>%
    mutate(n_samples = n()) %>%
    group_by(diet, NGS_pass) %>%
    summarize(rate = n() / first(n_samples) * 100, 
              n_samples = first(n_samples), 
              .groups='drop') %>%
    filter(NGS_pass == TRUE) %>%
    mutate(category = 'Diet') %>%
    rename('feature' = diet)
sr_diet
# overall success rate
sr_digsys = metadata %>%
    mutate(digestive_system = gsub('fermentation activity ', '', digestive_system)) %>%
    group_by(digestive_system) %>%
    mutate(n_samples = n()) %>%
    group_by(digestive_system, NGS_pass) %>%
    summarize(rate = n() / first(n_samples) * 100, 
              n_samples = first(n_samples), 
              .groups='drop') %>%
    filter(NGS_pass == TRUE) %>%
    mutate(category = 'Digestive system') %>%
    rename('feature' = digestive_system)
sr_digsys
# overall success rate
sr_hab = metadata %>%
    group_by(habitat) %>%
    mutate(n_samples = n()) %>%
    group_by(habitat, NGS_pass) %>%
    summarize(rate = n() / first(n_samples) * 100, 
              n_samples = first(n_samples), 
              .groups='drop') %>%
    filter(NGS_pass == TRUE) %>%
    mutate(category = 'Habitat') %>%
    rename('feature' = habitat)
sr_hab
# overall success rate
sr_cw = metadata %>%
    group_by(captive_wild) %>%
    mutate(n_samples = n()) %>%
    group_by(captive_wild, NGS_pass) %>%
    summarize(rate = n() / first(n_samples) * 100, 
              n_samples = first(n_samples), 
              .groups='drop') %>%
    filter(NGS_pass == TRUE) %>%
    mutate(category = 'Wild/captive') %>%
    rename('feature' = captive_wild)
sr_cw
sr = list(sr_class, sr_st, sr_diet, sr_digsys, sr_cw) %>%
    do.call(rbind, .)
df.dims(20)
sr
df.dims()
p = sr %>%
    ggplot(aes(feature, rate, fill=n_samples)) +
    geom_bar(stat='identity') +
    scale_fill_gradient('No. of\nsamples', low='blue', high='orange') +
    labs(x='Category', y='Success rate (%)') +
    facet_wrap(~ category, scales='free_x') +
    theme_bw() +
    theme(
        axis.text.x = element_text(angle=55, hjust=1)
    )

p.dims(7,5.5)
plot(p)
p = sr %>%
    filter(category %in% c('Class', 'Diet', 'Sample type', 'Wild/captive')) %>%
    ggplot(aes(feature, rate, fill=n_samples)) +
    geom_bar(stat='identity') +
    scale_fill_gradient('No. of\nsamples', low='blue', high='orange') +
    labs(x='Category', y='Success rate (%)') +
    facet_wrap(~ category, scales='free_x') +
    theme_bw() +
    theme(
        axis.text.x = element_text(angle=55, hjust=1)
    )

p.dims(6,5.5)
plot(p)
# overall success rate
sr_order = metadata %>%
    filter(class == 'Mammalia') %>%
    group_by(order) %>%
    mutate(n_samples = n()) %>%
    group_by(order, NGS_pass) %>%
    summarize(rate = n() / first(n_samples) * 100, 
              n_samples = first(n_samples), 
              .groups='drop') %>%
    filter(NGS_pass == TRUE) %>%
    mutate(category = 'Order') %>%
    rename('feature' = order)
sr_order
p = sr_order %>%
    mutate(feature = feature %>% reorder(-rate)) %>%
    ggplot(aes(feature, rate, fill=n_samples)) +
    geom_bar(stat='identity') +
    scale_fill_gradient('No. of\nsamples', low='blue', high='orange') +
    labs(x='Host order', y='Success rate (%)') +
    theme_bw() +
    theme(
        axis.text.x = element_text(angle=55, hjust=1)
    )

p.dims(5,3)
plot(p)<jupyter_output><empty_output><jupyter_text>Writing table of samples<jupyter_code>out_file = file.path(work_dir, 'NGS_pass-fail_metadata.tsv')
metadata %>%
    dplyr::select(SampleID, class, order, family, genus, scientific_name,
                  diet, sample_type, NGS_pass) %>%
    mutate(NGS_pass = ifelse(NGS_pass == TRUE, 'pass', 'fail')) %>%
    arrange(class, order, family, genus, scientific_name) %>%
    write_table(out_file)<jupyter_output>File written: /ebio/abt3_projects/Georg_animal_feces/data/16S_arch/NGS//NGS_pass-fail_metadata.tsv<jupyter_text>No. of samples<jupyter_code># number of sample samples
p = metadata %>%
    group_by(class, NGS_pass) %>%
    summarize(N = n(), .groups='drop') %>%
    ggplot(aes(class, N, fill=NGS_pass)) +
    geom_bar(stat='identity') +
    scale_fill_discrete('In final\nGS dataset?') +
    labs(x='Host class', y='No. of samples') +
    theme_bw() +
    theme(
        axis.text.x = element_text(angle=45, hjust=1)
    )

p.dims(4, 3)
p
# number of sample samples
p = metadata %>%
    group_by(sample_type, NGS_pass) %>%
    summarize(N = n(), .groups='drop') %>%
    ggplot(aes(sample_type, N, fill=NGS_pass)) +
    geom_bar(stat='identity') +
    scale_fill_discrete('In final\nGS dataset?') +
    labs(x='Sample type', y='No. of samples') +
    theme_bw() +
    theme(
        axis.text.x = element_text(angle=45, hjust=1)
    )

p.dims(4, 3)
p<jupyter_output><empty_output><jupyter_text>Passed per species<jupyter_code># samples per individual
x = metadata %>%
    filter(!is.na(scientific_name)) %>%
    group_by(scientific_name) %>%
    summarize(n_samples = n(), .groups='drop') 

# species with >1 sample
x %>% filter(n_samples > 1) %>% nrow
# all species
x %>% nrow
# samples per individual
metadata %>%
    filter(!is.na(scientific_name)) %>%
    group_by(scientific_name) %>%
    summarize(n_samples = n(), .groups='drop') %>%
    .$n_samples %>% summary_x
# % of individuals passed per species
meta_s = metadata %>%    
    filter(!is.na(scientific_name)) %>%
    group_by(class, scientific_name, NGS_pass) %>%
    summarize(N = n(), .groups='drop') %>%
    spread(NGS_pass, N, fill=0) %>%
    mutate(perc_pass = `TRUE` / (`TRUE` + `FALSE`) * 100)
meta_s$perc_pass %>% summary_x
meta_s
# plotting
p = meta_s %>%
    group_by(class) %>%
    mutate(n_species = n()) %>%
    ungroup() %>%
    ggplot(aes(class, perc_pass, fill=n_species)) +
    geom_boxplot() +
    scale_fill_gradient('No. of\nspecies', low='blue', high='orange') +
    labs(x='Host class', y='Per-species\nsuccess rate') +
    theme_bw() +
    theme(
        axis.text.x = element_text(angle=45, hjust=1)
    )

p.dims(4.5, 3)
p
# % of individuals passed per species (just those with >1 sample)
meta_s = metadata %>%    
    filter(!is.na(scientific_name)) %>%
    group_by(scientific_name) %>%
    mutate(n_samples = n()) %>%
    ungroup() %>%
    filter(n_samples > 1) %>%
    group_by(class, scientific_name, NGS_pass) %>%
    summarize(N = n(), .groups='drop') %>%
    spread(NGS_pass, N, fill=0) %>%
    mutate(perc_pass = `TRUE` / (`TRUE` + `FALSE`) * 100)
meta_s$perc_pass %>% summary_x
meta_s<jupyter_output><empty_output><jupyter_text>Mapping onto host tree<jupyter_code>overlap(meta_s$scientific_name, host_tree$tip.label)
overlap(meta_s$scientific_name, host_tree$tip.label, to_return='diff_x')
# filtering the tree
to_rm = setdiff(host_tree$tip.label, meta_s$scientific_name)
host_tree_f = ape::drop.tip(host_tree, to_rm)
host_tree_f
# creating itol mapping
df = meta_s %>%
    dplyr::select(scientific_name, `FALSE`, `TRUE`) %>%
    rename('passed' = `TRUE`,
           'failed' = `FALSE`) %>%
    as.data.frame
rownames(df) = df$scientific_name
df$scientific_name = NULL

F = file.path(work_dir, 'NGS_pass-fail_itol-bar.txt')
itol_multibar(df, 'NGS_pass_fail', F)
# writing out tree
F = file.path(work_dir, 'NGS_pass-fail_SpecD.nwk')
write.tree(host_tree_f, F)
cat('File written:', F, '\n')<jupyter_output>File written: /ebio/abt3_projects/Georg_animal_feces/data/16S_arch/NGS//NGS_pass-fail_SpecD.nwk<jupyter_text>Host metadata  sessionInfo<jupyter_code>sessionInfo()<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.15525749147676296, 'fraction_numerical': 0.06527005203660506, 'mean_word_length': 5.20456443083774, 'pattern_counts': {'":': 0, '<': 42, '<?xml version=': 0, '>': 192, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 2, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '11', 'dir': 'jupyter-structured-clean-dedup', 'id': '492749', 'license': 'MIT', 'n_tokens_mistral': 10449, 'n_tokens_neox': 8971, 'n_words': 1257, 'path': '00_misc/01_NGS_pass-fail.ipynb', 'repo_name': 'leylabmpi/16S-arc_vertebrate_paper', 'size': '22292'}