CAMISIM: Generate the abundance profiles

[ndreey]

With the chosen genomes, generate the abundance profiles. Smartest seems to create a bash script that generates it for us. Here is info from my comment on [#13]

HOW I UNDERSTAND THE ABUNDANCE CALCULATION The abundance is calculated based on the total sum of genome sizes.

• Say i have G1, G2, G3, ..., G10, Orchid genomes. And i want Orchid to have an abundance of 50%.
• G1:G5 is 1000bp each, G6:G10 is 1500bp each and Orchid is 12000bp

  Genome    Size
  G1        1000
  G2        1000
  ...   
  G6        1500
  G7        1500
  ...    
  Orchid    12000

• Calculate the total genome size
  • tot = 1000 x 5 + 1500 x 5 + 12000 = 19500bp
• Calculate the abundance value for each genome.
  • abu = 1 / (number of genomes - 1)
  • For G1 to G10 there are 10-1 = 9 genomes.
  • abu for G1 to G10 are 1/9 = 0.1111
• Set the abundance value of Orchid to 0.5
• Calculate the total abundance value for all genomes.
  • abu_tot = abu_orchid + sum(abu_G1:abu_G10)
  • abu_tot= 0.5 + 9 x 0.1111 = 1.5
• Normalize the abundance values so they sum up to 1.
  • nrm_abu_orchid = abu_ochid / abu_tot= 0.5 / 1.5 = 0.3333
  • nrm_abu_Gn = 0.1111 / 1.5 = 0.0741
• BOOOM there is your relative abundance. But NOTE there should not be a heading row in the abundance.tsv file
• However, it seems that abundance don't have to sum up to 1 as can be seen in the example above. But doing it this way i am able to sum all abundances to 1 "ish". 0.0741 x 10 + 0.3333 ~ 1

  Genome    Abundance
  G1        0.0741
  G2        0.0741
  ... 
  G10       0.0741
  Orchid    0.3333

Good info on these issues

• Can I use CAMISIM to simulate a sample with multiple strains of the same species?
• De novo community design: Genomes download
• ERROR: [MetagenomeSimulationPipeline] could not convert string to float: 'abundance' in line 80
• Guidance on abundance.tsv file and output
• custom abundance profile doesn't reflect abundance in output
• Generate simulated samples from a genome abundance table

[ndreey]

CREATING DATAFRAME

• I want a df with these columns genome_id size taxid tax_group group

pseudo code

list_genome_id = []
list_size = []
list_NCBI = []
list_tax_group = []
list_group = []

for *.fasta in source_genomes/:

    # Get genome_id
    match genome_id with *.fasta filename using genome_to_id.txt
    add genome_id to list_genome_id

    # Get size  
    match size with *.fasta filename using report_genome.txt
    add size to list_size

    # Get taxid
    match NCBI_ID with genome_id using metadata.tsv
    add NCBI_ID to list_NCBI

    # Get taxonomy group id
    match tax_group with NCBI_ID using taxonomic_profile.tsv
        in $TAXPATH get first number ^[0-9]|      # 2|1239|186801   --> 2
    add tax_group to list_tax_group

    # Get humanized group name
    match tax_group with group using if
    if tax_group != 2759:
        group = "not_euk"
    else:
        group = "euk"
    add group to list_group

    # Create dataframe
    df <- data.frame(genome_id = list_genome_id, size=, taxid=, tax_group=, group=)

[ndreey]

• Thaliana genome (GCF_000001735.4_TAIR10.1_genomic) has been replaced with Platanthera_zijinensis_chr

• Added these to source_genomes/

  • Ceratobasidium_sp_CerAGI
  • Rhizoctonia_solani_Rhisola1
  • Tulasnella_calospora_Tulcal1

[ndreey]

26

[ndreey]

Going with 70-80% Fungi, 20-30% Bacteria/Archaea and 0.5-2% Plasmids/Circular DNA/Virus #17

[ndreey]

Example

• Lets generate 10 GB data from the total genome size of all the 100 genomes:
• Host abundance: 50% --> 5GB data belongs to host
• Endophyte abundance: 1 - Host abundance --> 5GB belong to endophytes.

The relative abundance between the endophytes are:

• Fungi: 67,8% | 3.39GB
  • OMF: 45,2% | 1.53GB (2 x rfungi)
    • Three OMFs --> 45.2%/3 = 15.06667% | 0.2305GB
  • rfungi: 22.6% |
    • Decided randomly but gotta sum up to 22.6%
• Bark: 29,7%
  • Decided randomly but gotta sum up to 29.7%
• Plasm: 2.5%
  • Decided randomly but gotta sum up to 2.5%