Goal:

  • Download most up-to-date version of NCBI 'complete' genomes

Setting variables


In [5]:
workDir = '/var/seq_data/ncbi_db/genome/Jan2016/'

Init


In [6]:
import os
%load_ext rpy2.ipython

In [60]:
%%R
library(ggplot2)
library(dplyr)
library(tidyr)
library(genomes)
library(permute)

In [8]:
if not os.path.isdir(workDir):
    os.makedirs(workDir)

Assessing NCBI prokaryote genome listing


In [11]:
%%R
data(proks)
summary(proks)


$`Total genomes`
[1] 27570 genome projects on Sep 04, 2014

$`By status`
                     Total
Contig               13074
Scaffold             10718
Gapless Chromosome    3053
Chromosome             373
Chromosome with gaps   343
Complete                 9

$`Recent submissions`
  released   name                      status  
1 2014-09-02 Altuibacter lentus        Scaffold
2 2014-09-02 Bacillus cereus ATCC 4342 Scaffold
3 2014-09-02 Bacillus licheniformis    Scaffold
4 2014-09-02 Bacillus megaterium       Scaffold
5 2014-09-02 Paenibacillus macerans    Scaffold


In [12]:
%%R
update(proks)
summary(proks)


/opt/anaconda/lib/python2.7/site-packages/rpy2/robjects/functions.py:106: UserWarning: proks has been successfully updated

  res = super(Function, self).__call__(*new_args, **new_kwargs)
/opt/anaconda/lib/python2.7/site-packages/rpy2/robjects/functions.py:106: UserWarning: 30013 new project IDs added

  res = super(Function, self).__call__(*new_args, **new_kwargs)
/opt/anaconda/lib/python2.7/site-packages/rpy2/robjects/functions.py:106: UserWarning: 237 old project IDs removed

  res = super(Function, self).__call__(*new_args, **new_kwargs)
$`Total genomes`
[1] 57892 genome projects on Jan 12, 2016

$`By status`
                Total
Contig          35209
Scaffold        17080
Complete Genome  4732
Chromosome        871

$`Recent submissions`
  released   name                 status         
1 2016-01-07 Campylobacter jejuni Complete Genome
2 2016-01-07 Campylobacter jejuni Complete Genome
3 2016-01-07 Campylobacter jejuni Complete Genome
4 2016-01-07 Campylobacter jejuni Complete Genome
5 2016-01-07 Campylobacter jejuni Complete Genome


In [13]:
%%R -w 600 -h 300
# plotting GC distribution
ggplot(proks, aes(gc)) +
    geom_histogram(binwidth=1) +
    geom_vline(xintercept=50, linetype='dashed', color='red', alpha=0.7) +
    labs(x='G+C') +
    theme(
        text = element_text(size=16)
    )


Complete genomes


In [14]:
%%R
proks.complete = proks %>% as.data.frame %>% 
    filter(status == 'Complete Genome') 
proks.complete %>% head(n=3)


     pid                           name          status   released  taxid
1  12997 Acaryochloris marina MBIC11017 Complete Genome 2007-10-16 329726
2  60713 Acetobacterium woodii DSM 1030 Complete Genome 2012-02-14 931626
3 242487       Acetobacter pasteurianus Complete Genome 2015-07-21    438
   bioproject          group              subgroup    size      gc
1  PRJNA12997  Cyanobacteria Oscillatoriophycideae 8.36160 46.9889
2  PRJNA60713     Firmicutes            Clostridia 4.04478 39.3000
3 PRJNA242487 Proteobacteria   Alphaproteobacteria 2.80615 53.3000
         refseq      insdc
1   NC_009925.1 CP000828.1
2   NC_016894.1 CP002987.1
3 NZ_CP012111.1 CP012111.1
                                                                                               plasmid.refseq
1 NC_009930.1,NC_009928.1,NC_009927.1,NC_009931.1,NC_009933.1,NC_009934.1,NC_009932.1,NC_009926.1,NC_009929.1
2                                                                                                        <NA>
3                                                                                                        <NA>
                                                                                       plasmid.insdc
1 CP000842.1,CP000840.1,CP000839.1,CP000843.1,CP000845.1,CP000846.1,CP000844.1,CP000838.1,CP000841.1
2                                                                                               <NA>
3                                                                                               <NA>
   wgs scaffolds genes proteins   modified                             center
1 <NA>        10  7469     7187 2015-07-30              Washington University
2 <NA>         1  3649     3521 2015-08-18 Georg-August-University Goettingen
3 <NA>         1  2688     2535 2015-12-09      Zhejiang Gongshang University
     biosample        assembly reference ftp   pubmed
1 SAMN02604308 GCA_000018105.1      REFR  NA 18252824
2 SAMN02603267 GCA_000247605.1      REFR  NA 22479398
3 SAMN02709032 GCA_001183745.1      <NA>  NA     <NA>

In [15]:
%%R -w 600 -h 300
# plotting GC distribution
ggplot(proks.complete, aes(gc)) +
    geom_histogram(binwidth=1) +
    geom_vline(xintercept=50, linetype='dashed', color='red', alpha=0.7) +
    labs(x='G+C') +
    theme(
        text = element_text(size=16)
    )



In [16]:
%%R
proks.complete$gc %>% summary


   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  13.50   37.67   46.70   47.77   58.80   74.90 

Adding taxonomy

  • using taxID

In [17]:
prokFile = os.path.join(workDir, 'proks_complete.txt')

In [18]:
%%R -i prokFile

write.table(proks.complete, prokFile, sep='\t', row.names=F, quote=F)
cat('File written: ', prokFile, '\n')


File written:  /var/seq_data/ncbi_db/genome/Jan2016/proks_complete.txt 

In [20]:
taxFile = os.path.splitext(prokFile)[0] + '_tax.txt'

In [156]:
!cd $workDir; \
    tail -n +2 $prokFile | \
    cut -f 5 | sort -u | \
    seqDB_tools taxid2taxonomy -p 30 > $taxFile


sh: 0: getcwd() failed: No such file or directory
Failed query (try=1): 738
Failed query (try=1): 741091
Failed query (try=1): 741093
Failed query (try=1): 91891

Reading in files


In [21]:
%%R -i taxFile

df.tax = read.delim(taxFile, sep='\t') %>%
    distinct(taxid)
df.proks.complete = dplyr::inner_join(proks.complete, df.tax, c('taxid' = 'taxid'))

# checking join
proks.complete %>% nrow %>% print
df.proks.complete %>% nrow %>% print
df.proks.complete %>% head(n=3)


[1] 4732
[1] 4732
     pid                           name          status   released  taxid
1  12997 Acaryochloris marina MBIC11017 Complete Genome 2007-10-16 329726
2  60713 Acetobacterium woodii DSM 1030 Complete Genome 2012-02-14 931626
3 242487       Acetobacter pasteurianus Complete Genome 2015-07-21    438
   bioproject          group              subgroup    size      gc
1  PRJNA12997  Cyanobacteria Oscillatoriophycideae 8.36160 46.9889
2  PRJNA60713     Firmicutes            Clostridia 4.04478 39.3000
3 PRJNA242487 Proteobacteria   Alphaproteobacteria 2.80615 53.3000
         refseq      insdc
1   NC_009925.1 CP000828.1
2   NC_016894.1 CP002987.1
3 NZ_CP012111.1 CP012111.1
                                                                                               plasmid.refseq
1 NC_009930.1,NC_009928.1,NC_009927.1,NC_009931.1,NC_009933.1,NC_009934.1,NC_009932.1,NC_009926.1,NC_009929.1
2                                                                                                        <NA>
3                                                                                                        <NA>
                                                                                       plasmid.insdc
1 CP000842.1,CP000840.1,CP000839.1,CP000843.1,CP000845.1,CP000846.1,CP000844.1,CP000838.1,CP000841.1
2                                                                                               <NA>
3                                                                                               <NA>
   wgs scaffolds genes proteins   modified                             center
1 <NA>        10  7469     7187 2015-07-30              Washington University
2 <NA>         1  3649     3521 2015-08-18 Georg-August-University Goettingen
3 <NA>         1  2688     2535 2015-12-09      Zhejiang Gongshang University
     biosample        assembly reference ftp   pubmed superkingdom
1 SAMN02604308 GCA_000018105.1      REFR  NA 18252824     Bacteria
2 SAMN02603267 GCA_000247605.1      REFR  NA 22479398     Bacteria
3 SAMN02709032 GCA_001183745.1      <NA>  NA     <NA>     Bacteria
          phylum               class            order           family
1  Cyanobacteria                <NA>    Chroococcales             <NA>
2     Firmicutes          Clostridia    Clostridiales   Eubacteriaceae
3 Proteobacteria Alphaproteobacteria Rhodospirillales Acetobacteraceae
           genus               species
1  Acaryochloris  Acaryochloris marina
2 Acetobacterium Acetobacterium woodii
3    Acetobacter                  <NA>

Just Bacteria


In [217]:
%%R
df.bac.complete = df.proks.complete %>%
    filter(superkingdom == 'Bacteria')

df.bac.complete %>% nrow


[1] 4498

Phylum representation


In [218]:
%%R -w 800
df.bac.complete.s = df.bac.complete %>%
    group_by(phylum) %>%
    summarize(n = n())

ggplot(df.bac.complete.s, aes(phylum, n)) +
    geom_bar(stat='identity') +
    scale_y_log10() +
    labs(y = 'Number of genomes') +
    theme_bw() +
    theme(
        text = element_text(size=16),
        axis.text.x = element_text(angle=60, hjust=1)
    )


Other Filtering

Genome size


In [219]:
%%R
df.bac.complete %>%
    filter(size < 0.8) %>%
    arrange(size) %>%
    select(size, name) %>% tail(n=20)


        size                                                         name
102 0.723970                Aster yellows witches'-broom phytoplasma AYWB
103 0.727289                                  Ureaplasma parvum serovar 3
104 0.742431                                Mycoplasma suis str. Illinois
105 0.749321  Blochmannia endosymbiont of Polyrhachis (Hedomyrma) turneri
106 0.751679                  Ureaplasma parvum serovar 3 str. ATCC 27815
107 0.751719                 Ureaplasma parvum serovar 3 str. ATCC 700970
108 0.752630  Parcubacteria (Campbellbacteria) bacterium GW2011_OD1_34_28
109 0.754729                     Candidatus Ishikawaella capsulata Mpkobe
110 0.756845                 Candidatus Mycoplasma haemolamae str. Purdue
111 0.759425                         Candidatus Baumannia cicadellinicola
112 0.760240                                  Mycoplasma bovoculi M165/69
113 0.773940 Blochmannia endosymbiont of Camponotus (Colobopsis) obliquus
114 0.777079                                       Mycoplasma mobile 163K
115 0.778866                             Mycoplasma flocculare ATCC 27399
116 0.791219                  Candidatus Blochmannia chromaiodes str. 640
117 0.791654              Candidatus Blochmannia pennsylvanicus str. BPEN
118 0.793224                                         Mesoplasma florum L1
119 0.793841                                      Mycoplasma californicum
120 0.799088                             Mycoplasma californicum HAZ160_1
121 0.799476                                       Mycoplasma synoviae 53

In [220]:
%%R

cat('Pre-filter:', df.bac.complete %>% nrow, '\n')
df.bac.complete = df.bac.complete %>%
    filter(size > 0.8) 
cat('Post-filter:', df.bac.complete %>% nrow, '\n')


Pre-filter: 4498 
Post-filter: 4377 

removing what are really phage/plasmid genomes


In [221]:
%%R
cat('Pre-filter:', df.bac.complete %>% nrow, '\n')

to.rm = c("Thermoanaerobacterium saccharolyticum JW/SL-YS485",
          "Streptococcus salivarius 57.I")

df.bac.complete = df.bac.complete %>%
    filter(! name %in% to.rm)

cat('Post-filter:', df.bac.complete %>% nrow, '\n')


Pre-filter: 4377 
Post-filter: 4375 

Random selection of just 1 taxon per species


In [222]:
%%R

# parsing out representatives
df.bac.complete.rep = df.bac.complete %>%
    filter(! is.na(refseq)) %>%
    group_by(species) %>%
    sample_n(1) %>%
    ungroup()

cat('Number of representative genomes:', df.bac.complete.rep %>% nrow, '\n')


Number of representative genomes: 1234 

Phylum representation


In [223]:
%%R -w 800
df.bac.complete.rep.s = df.bac.complete.rep %>%
    group_by(phylum) %>%
    summarize(n = n())

ggplot(df.bac.complete.rep.s, aes(phylum, n)) +
    geom_bar(stat='identity') +
    labs(y = 'Number of representative genomes') +
    theme_bw() +
    theme(
        text = element_text(size=16),
        axis.text.x = element_text(angle=60, hjust=1)
    )


genome GC content


In [224]:
%%R -w 600 -h 300
# plotting GC distribution
ggplot(df.bac.complete.rep, aes(gc)) +
    geom_histogram(binwidth=1) +
    geom_vline(xintercept=50, linetype='dashed', color='red', alpha=0.7) +
    labs(x='G+C') +
    theme(
        text = element_text(size=16)
    )


Notes

  • A very bi-modal distribution (as see in Youngblut & Buckley (2014))

Reducing bias toward low G+C endosymbionts


In [225]:
%%R -h 400
# culmulative function of Genomes w/ GC

df.bac.complete.rep.f = df.bac.complete.rep %>%
    group_by() %>%
    mutate(total_genomes = n()) %>%
    ungroup() %>%
    arrange(gc) %>%
    mutate(n = 1,
           culm_perc_total = cumsum(n) / total_genomes * 100) %>%
    select(gc, culm_perc_total, phylum, genus) 
    

# plot
ggplot(df.bac.complete.rep.f, aes(gc, culm_perc_total)) +
    geom_point() +
    geom_line() +
    labs(x='Genome G+C', y='Culmulative total (%)') +
    theme_bw() +
    theme(
        text = element_text(size=16)
    )



In [226]:
%%R -w 1100 -h 500
GC_cutoff = 30

df.bac.complete.rep.f = df.bac.complete.rep %>% 
    group_by(genus) %>%
    mutate(genus_medianGC = median(gc)) %>%
    ungroup() %>%
    filter(genus_medianGC <= GC_cutoff) %>%
    select(phylum, genus, gc)

df.bac.complete.rep.f$genus = reorder(df.bac.complete.rep.f$genus, df.bac.complete.rep.f$phylum)

ggplot(df.bac.complete.rep.f, aes(genus, gc, color=phylum)) +
    geom_boxplot() +
    theme_bw() +
    theme(
        text = element_text(size=16),
        axis.text.x = element_text(angle=60, hjust=1)
    )


Making low GC genomes just 1% of dataset


In [228]:
%%R
df.bac.complete.rep = df.bac.complete.rep %>% 
    group_by(genus) %>%
    mutate(genus_medianGC = median(gc)) %>%
    ungroup()

df.lowGC = df.bac.complete.rep %>%
    filter(genus_medianGC <= 30)

df.highGC = df.bac.complete.rep %>%
    filter(genus_medianGC > 30)

n.lowGC = df.lowGC %>% nrow
n.total = df.bac.complete.rep %>% nrow

cat('Number of low GC genomes: ', n.lowGC, '\n')
cat('Low GC genomes make up ', n.lowGC / n.total * 100, '% of the dataset\n')


Number of low GC genomes:  79 
Low GC genomes make up  6.401945 % of the dataset

In [229]:
%%R -h 300 -w 600
n.sample = round(n.total * 0.01, 0)
cat('Sampling', n.sample, 'low GC genomes\n')

to.keep = sample(df.bac.complete.rep$name, n.sample, replace=FALSE)
to.keep = append(to.keep, df.highGC$name)

df.bac.complete.rep.p = df.bac.complete.rep %>%
    filter(name %in% to.keep)    

cat('Number of representative genomes:', df.bac.complete.rep.p %>% nrow, '\n')

# plotting GC distribution
ggplot(df.bac.complete.rep.p, aes(gc)) +
    geom_histogram(binwidth=1) +
    geom_vline(xintercept=50, linetype='dashed', color='red', alpha=0.7) +
    labs(x='G+C') +
    theme(
        text = element_text(size=16)
    )


Sampling 12 low GC genomes
Number of representative genomes: 1155 

Sequence download


In [230]:
%%R -i workDir

outFile = file.path(workDir, 'bac_complete_spec-rep1.txt')
write.table(df.bac.complete.rep.p, outFile, sep='\t', quote=F, row.names=F)

In [ ]:
!cd $workDir; \
    seqDB_tools accession-GI2fasta \
        -a 11 -n 2 -f 30 -header \
        -o bac_complete_spec-rep1 \
        < bac_complete_spec-rep1.txt \
        2> bac_complete_spec-rep1.log

In [239]:
!cd $workDir; \
    seqDB_tools accession-GI2fasta \
        -a 11 -n 2 -f 30 -b 1 -header \
        -o bac_complete_spec-rep1_tmp \
        < re_bac.txt


Starting batch->trial: 0->1
Starting batch->trial: 1->1
Starting batch->trial: 3->1
Starting batch->trial: 2->1

--------------------- WARNING ---------------------
MSG: No whitespace allowed in FASTA ID [NC_015387|Marinithermus hydrothermalis DSM 14884, complete genome.]
---------------------------------------------------

--------------------- WARNING ---------------------
MSG: No whitespace allowed in FASTA ID [NC_015387|Marinithermus hydrothermalis DSM 14884, complete genome.]
---------------------------------------------------

--------------------- WARNING ---------------------
MSG: No whitespace allowed in FASTA ID [NC_017094|Leptospirillum ferrooxidans C2-3 DNA, complete genome.]
---------------------------------------------------

--------------------- WARNING ---------------------
MSG: No whitespace allowed in FASTA ID [NC_017094|Leptospirillum ferrooxidans C2-3 DNA, complete genome.]
---------------------------------------------------

--------------------- WARNING ---------------------
MSG: No whitespace allowed in FASTA ID [NC_022041|Paracoccus aminophilus JCM 7686, complete genome.]
---------------------------------------------------

--------------------- WARNING ---------------------
MSG: No whitespace allowed in FASTA ID [NC_022041|Paracoccus aminophilus JCM 7686, complete genome.]
---------------------------------------------------

--------------------- WARNING ---------------------
MSG: No whitespace allowed in FASTA ID [NC_020908|Octadecabacter arcticus 238, complete genome.]
---------------------------------------------------

--------------------- WARNING ---------------------
MSG: No whitespace allowed in FASTA ID [NC_020908|Octadecabacter arcticus 238, complete genome.]
---------------------------------------------------

In [ ]:
Leptospirillum_ferrooxidans_C2-3.fna
bac_complete_spec-rep1/Marinithermus_hydrothermalis_DSM_14884.fna
bac_complete_spec-rep1/Octadecabacter_arcticus_238.fna
bac_complete_spec-rep1/Paracoccus_aminophilus_JCM_7686.fna

Checking output


In [240]:
genomeDir = os.path.join(workDir, 'bac_complete_spec-rep1')
# number of genomes downloaded
!printf "Number of bacterial genomes: "
!cd $genomeDir; \
    find . -name "*.fna" | wc -l


Number of bacterial genomes: 1147

In [243]:
# file size
!echo "Genome file size distribution (bytes):"
!cd $genomeDir; \
    ls -tlc *.fna | \
    perl -pe 's/ +/\t/g' | \
    cut -f 5 | NY_misc_perl stats_descriptive


Genome file size distribution (bytes):
1	min	721810.00
1	Q1	2386016.00
1	mean	3914102.84
1	median	3668349.00
1	Q3	4991326.00
1	max	15028601.00
1	stdev	1945831.55

In [242]:
# checking for non-bacterial genomes
!find $genomeDir -name "*fna" | xargs -P 20 egrep "phage|virus|phage"

Renaming genomes


In [251]:
genomeDirRn = genomeDir + '_rn'
genomeDirRn


Out[251]:
'/var/seq_data/ncbi_db/genome/Jan2016/bac_complete_spec-rep1_rn'

In [252]:
# making sure each sequence is unique
!cd $genomeDir; \
    find . -name "*fna" | \
    SIPSim genome_rename -n 24 --prefix $genomeDirRn -


File written: /var/seq_data/ncbi_db/genome/Jan2016/bac_complete_spec-rep1_rn/Chlamydia_pecorum_P787.fna
File written: /var/seq_data/ncbi_db/genome/Jan2016/bac_complete_spec-rep1_rn/Streptococcus_pyogenes_MGAS8232.fna
File written: /var/seq_data/ncbi_db/genome/Jan2016/bac_complete_spec-rep1_rn/Rickettsia_australis_str_Cutlack.fna
File written: /var/seq_data/ncbi_db/genome/Jan2016/bac_complete_spec-rep1_rn/Chlamydia_avium_10DC88.fna
File written: /var/seq_data/ncbi_db/genome/Jan2016/bac_complete_spec-rep1_rn/Haemophilus_ducreyi_35000HP.fna
File written: /var/seq_data/ncbi_db/genome/Jan2016/bac_complete_spec-rep1_rn/Thermotoga_maritima_MSB8.fna
File written: /var/seq_data/ncbi_db/genome/Jan2016/bac_complete_spec-rep1_rn/Lactococcus_piscium_MKFS47.fna
File written: /var/seq_data/ncbi_db/genome/Jan2016/bac_complete_spec-rep1_rn/Corynebacterium_falsenii_DSM_44353.fna
File written: /var/seq_data/ncbi_db/genome/Jan2016/bac_complete_spec-rep1_rn/Helicobacter_cinaedi_PAGU611.fna
File written: /var/seq_data/ncbi_db/genome/Jan2016/bac_complete_spec-rep1_rn/Thermodesulfobacterium_commune_DSM_2178.fna
File written: /var/seq_data/ncbi_db/genome/Jan2016/bac_complete_spec-rep1_rn/Desulfotomaculum_reducens_MI-1.fna
File written: /var/seq_data/ncbi_db/genome/Jan2016/bac_complete_spec-rep1_rn/Bacillus_megaterium_WSH-002.fna
File written: /var/seq_data/ncbi_db/genome/Jan2016/bac_complete_spec-rep1_rn/Pantoea_ananatis_AJ13355.fna
File written: /var/seq_data/ncbi_db/genome/Jan2016/bac_complete_spec-rep1_rn/Saccharophagus_degradans_2-40.fna
File written: /var/seq_data/ncbi_db/genome/Jan2016/bac_complete_spec-rep1_rn/Simkania_negevensis_Z.fna
File written: /var/seq_data/ncbi_db/genome/Jan2016/bac_complete_spec-rep1_rn/Solibacillus_silvestris_StLB046.fna
File written: /var/seq_data/ncbi_db/genome/Jan2016/bac_complete_spec-rep1_rn/Bdellovibrio_bacteriovorus_str_Tiberius.fna
File written: /var/seq_data/ncbi_db/genome/Jan2016/bac_complete_spec-rep1_rn/Streptococcus_uberis_0140J.fna
File written: /var/seq_data/ncbi_db/genome/Jan2016/bac_complete_spec-rep1_rn/Pseudomonas_mendocina_ymp.fna
File written: /var/seq_data/ncbi_db/genome/Jan2016/bac_complete_spec-rep1_rn/Rickettsia_bellii_OSU_85-389.fna
File written: /var/seq_data/ncbi_db/genome/Jan2016/bac_complete_spec-rep1_rn/Desulfovibrio_magneticus_RS-1.fna
File written: /var/seq_data/ncbi_db/genome/Jan2016/bac_complete_spec-rep1_rn/Segniliparus_rotundus_DSM_44985.fna
File written: /var/seq_data/ncbi_db/genome/Jan2016/bac_complete_spec-rep1_rn/Sphaerobacter_thermophilus_DSM_20745.fna
File written: /var/seq_data/ncbi_db/genome/Jan2016/bac_complete_spec-rep1_rn/Thioalkalivibrio_sulfidiphilus_HL-EbGr7.fna
File written: /var/seq_data/ncbi_db/genome/Jan2016/bac_complete_spec-rep1_rn/Pseudomonas_pseudoalcaligenes_CECT_5344.fna
File written: /var/seq_data/ncbi_db/genome/Jan2016/bac_complete_spec-rep1_rn/Rhodoferax_ferrireducens_T118.fna
File written: /var/seq_data/ncbi_db/genome/Jan2016/bac_complete_spec-rep1_rn/Bifidobacterium_kashiwanohense_JCM_15439_DSM_21854.fna
File written: /var/seq_data/ncbi_db/genome/Jan2016/bac_complete_spec-rep1_rn/Enterococcus_casseliflavus_EC20.fna
File written: /var/seq_data/ncbi_db/genome/Jan2016/bac_complete_spec-rep1_rn/Thioflavicoccus_mobilis_8321.fna
File written: /var/seq_data/ncbi_db/genome/Jan2016/bac_complete_spec-rep1_rn/Eubacterium_eligens_ATCC_27750.fna
File written: /var/seq_data/ncbi_db/genome/Jan2016/bac_complete_spec-rep1_rn/Providencia_stuartii_MRSN_2154.fna
File written: /var/seq_data/ncbi_db/genome/Jan2016/bac_complete_spec-rep1_rn/Sphingomonas_wittichii_RW1.fna
File written: /var/seq_data/ncbi_db/genome/Jan2016/bac_complete_spec-rep1_rn/Paludibacter_propionicigenes_WB4.fna
File written: /var/seq_data/ncbi_db/genome/Jan2016/bac_complete_spec-rep1_rn/Candidatus_Kinetoplastibacterium_galatii_TCC219.fna
File written: /var/seq_data/ncbi_db/genome/Jan2016/bac_complete_spec-rep1_rn/Mycobacterium_gilvum_Spyr1.fna
File written: /var/seq_data/ncbi_db/genome/Jan2016/bac_complete_spec-rep1_rn/Aeromonas_media_WS.fna
File written: /var/seq_data/ncbi_db/genome/Jan2016/bac_complete_spec-rep1_rn/Xanthomonas_oryzae_pv_oryzicola.fna
File written: /var/seq_data/ncbi_db/genome/Jan2016/bac_complete_spec-rep1_rn/Lactobacillus_buchneri_CD034.fna
File written: /var/seq_data/ncbi_db/genome/Jan2016/bac_complete_spec-rep1_rn/Streptococcus_oligofermentans_AS_1_3089.fna
File written: /var/seq_data/ncbi_db/genome/Jan2016/bac_complete_spec-rep1_rn/Thermoanaerobacter_mathranii_subsp_mathranii_str_A3.fna
File written: /var/seq_data/ncbi_db/genome/Jan2016/bac_complete_spec-rep1_rn/Rhodobacter_capsulatus_SB_1003.fna
File written: /var/seq_data/ncbi_db/genome/Jan2016/bac_complete_spec-rep1_rn/Francisella_noatunensis_subsp_orientalis_str_Toba_04.fna
File written: /var/seq_data/ncbi_db/genome/Jan2016/bac_complete_spec-rep1_rn/Acidothermus_cellulolyticus_11B.fna
File written: /var/seq_data/ncbi_db/genome/Jan2016/bac_complete_spec-rep1_rn/Acidaminococcus_fermentans_DSM_20731.fna
File written: /var/seq_data/ncbi_db/genome/Jan2016/bac_complete_spec-rep1_rn/Xenorhabdus_poinarii_G6.fna
File written: /var/seq_data/ncbi_db/genome/Jan2016/bac_complete_spec-rep1_rn/Mycobacterium_bovis_BCG.fna
File written: /var/seq_data/ncbi_db/genome/Jan2016/bac_complete_spec-rep1_rn/Lactobacillus_salivarius_UCC118.fna
File written: /var/seq_data/ncbi_db/genome/Jan2016/bac_complete_spec-rep1_rn/Leuconostoc_kimchii_IMSNU_11154.fna
File written: /var/seq_data/ncbi_db/genome/Jan2016/bac_complete_spec-rep1_rn/Candidatus_Kinetoplastibacterium_crithidii_ex_Angomonas_deanei_ATCC_30255.fna
File written: /var/seq_data/ncbi_db/genome/Jan2016/bac_complete_spec-rep1_rn/Syntrophobotulus_glycolicus_DSM_8271.fna
File written: /var/seq_data/ncbi_db/genome/Jan2016/bac_complete_spec-rep1_rn/Campylobacter_hominis_ATCC_BAA-381.fna
File written: /var/seq_data/ncbi_db/genome/Jan2016/bac_complete_spec-rep1_rn/Desulfovibrio_salexigens_DSM_2638.fna
File written: /var/seq_data/ncbi_db/genome/Jan2016/bac_complete_spec-rep1_rn/Acidiphilium_cryptum_JF-5.fna
File written: /var/seq_data/ncbi_db/genome/Jan2016/bac_complete_spec-rep1_rn/Pseudomonas_denitrificans_ATCC_13867.fna
File written: /var/seq_data/ncbi_db/genome/Jan2016/bac_complete_spec-rep1_rn/Kangiella_koreensis_DSM_16069.fna
File written: /var/seq_data/ncbi_db/genome/Jan2016/bac_complete_spec-rep1_rn/Adlercreutzia_equolifaciens_DSM_19450.fna
File written: /var/seq_data/ncbi_db/genome/Jan2016/bac_complete_spec-rep1_rn/Taylorella_equigenitalis_MCE9.fna
File written: /var/seq_data/ncbi_db/genome/Jan2016/bac_complete_spec-rep1_rn/Geobacillus_thermoleovorans_CCB_US3_UF5.fna
File written: /var/seq_data/ncbi_db/genome/Jan2016/bac_complete_spec-rep1_rn/Parvibaculum_lavamentivorans_DS-1.fna
File written: /var/seq_data/ncbi_db/genome/Jan2016/bac_complete_spec-rep1_rn/Acidithiobacillus_ferrooxidans_ATCC_23270.fna
File written: /var/seq_data/ncbi_db/genome/Jan2016/bac_complete_spec-rep1_rn/Stenotrophomonas_maltophilia_K279a.fna
File written: /var/seq_data/ncbi_db/genome/Jan2016/bac_complete_spec-rep1_rn/Thermobaculum_terrenum_ATCC_BAA-798.fna
File written: /var/seq_data/ncbi_db/genome/Jan2016/bac_complete_spec-rep1_rn/Vibrio_parahaemolyticus_RIMD_2210633.fna
File written: /var/seq_data/ncbi_db/genome/Jan2016/bac_complete_spec-rep1_rn/Leptospirillum_ferriphilum_ML-04.fna
File written: /var/seq_data/ncbi_db/genome/Jan2016/bac_complete_spec-rep1_rn/Streptomyces_collinus_Tu_365.fna
File written: /var/seq_data/ncbi_db/genome/Jan2016/bac_complete_spec-rep1_rn/Amycolicicoccus_subflavus_DQS3-9A1.fna
File written: /var/seq_data/ncbi_db/genome/Jan2016/bac_complete_spec-rep1_rn/Prevotella_sp_oral_taxon_299_str_F0039.fna
File written: /var/seq_data/ncbi_db/genome/Jan2016/bac_complete_spec-rep1_rn/Microbacterium_testaceum_StLB037.fna
File written: /var/seq_data/ncbi_db/genome/Jan2016/bac_complete_spec-rep1_rn/Mesotoga_prima_MesG1_Ag_4_2.fna
File written: /var/seq_data/ncbi_db/genome/Jan2016/bac_complete_spec-rep1_rn/Helicobacter_mustelae_12198.fna
File written: /var/seq_data/ncbi_db/genome/Jan2016/bac_complete_spec-rep1_rn/Sanguibacter_keddieii_DSM_10542.fna
File written: /var/seq_data/ncbi_db/genome/Jan2016/bac_complete_spec-rep1_rn/Corynebacterium_kroppenstedtii_DSM_44385.fna
File written: /var/seq_data/ncbi_db/genome/Jan2016/bac_complete_spec-rep1_rn/Bacteroides_vulgatus_ATCC_8482.fna
File written: /var/seq_data/ncbi_db/genome/Jan2016/bac_complete_spec-rep1_rn/Bifidobacterium_angulatum_DSM_20098_JCM_7096.fna
File written: /var/seq_data/ncbi_db/genome/Jan2016/bac_complete_spec-rep1_rn/Treponema_paraluiscuniculi_Cuniculi_A.fna
File written: /var/seq_data/ncbi_db/genome/Jan2016/bac_complete_spec-rep1_rn/Pseudomonas_balearica_DSM_6083.fna
File written: /var/seq_data/ncbi_db/genome/Jan2016/bac_complete_spec-rep1_rn/Desulfurispirillum_indicum_S5.fna
File written: /var/seq_data/ncbi_db/genome/Jan2016/bac_complete_spec-rep1_rn/Intrasporangium_calvum_DSM_43043.fna
File written: /var/seq_data/ncbi_db/genome/Jan2016/bac_complete_spec-rep1_rn/Lactococcus_garvieae_Lg2.fna
File written: /var/seq_data/ncbi_db/genome/Jan2016/bac_complete_spec-rep1_rn/Kosakonia_sacchari_SP1.fna
File written: /var/seq_data/ncbi_db/genome/Jan2016/bac_complete_spec-rep1_rn/Robiginitalea_biformata_HTCC2501.fna
File written: /var/seq_data/ncbi_db/genome/Jan2016/bac_complete_spec-rep1_rn/Nitrobacter_hamburgensis_X14.fna
File written: /var/seq_data/ncbi_db/genome/Jan2016/bac_complete_spec-rep1_rn/Streptococcus_constellatus_subsp_pharyngis_C1050.fna
File written: /var/seq_data/ncbi_db/genome/Jan2016/bac_complete_spec-rep1_rn/Pseudomonas_resinovorans_NBRC_106553.fna
File written: /var/seq_data/ncbi_db/genome/Jan2016/bac_complete_spec-rep1_rn/Cellulomonas_flavigena_DSM_20109.fna
File written: /var/seq_data/ncbi_db/genome/Jan2016/bac_complete_spec-rep1_rn/Roseburia_hominis_A2-183.fna
File written: /var/seq_data/ncbi_db/genome/Jan2016/bac_complete_spec-rep1_rn/Bacillus_thuringiensis_serovar_kurstaki_str_HD73.fna
File written: /var/seq_data/ncbi_db/genome/Jan2016/bac_complete_spec-rep1_rn/Anaplasma_marginale_str_St_Maries.fna
File written: /var/seq_data/ncbi_db/genome/Jan2016/bac_complete_spec-rep1_rn/Streptococcus_infantarius_subsp_infantarius_CJ18.fna
File written: /var/seq_data/ncbi_db/genome/Jan2016/bac_complete_spec-rep1_rn/Burkholderia_phytofirmans_PsJN.fna
File written: /var/seq_data/ncbi_db/genome/Jan2016/bac_complete_spec-rep1_rn/Fibrobacter_succinogenes_subsp_succinogenes_S85.fna
File written: /var/seq_data/ncbi_db/genome/Jan2016/bac_complete_spec-rep1_rn/Pectobacterium_wasabiae_WPP163.fna
File written: /var/seq_data/ncbi_db/genome/Jan2016/bac_complete_spec-rep1_rn/Brucella_ovis_ATCC_25840.fna
File written: /var/seq_data/ncbi_db/genome/Jan2016/bac_complete_spec-rep1_rn/Tropheryma_whipplei_str_Twist.fna
File written: /var/seq_data/ncbi_db/genome/Jan2016/bac_complete_spec-rep1_rn/Methylacidiphilum_fumariolicum_SolV.fna
**OUTPUT MUTED**

Checking genome info

  • Making sure there are no outliers that need to be filtered out

In [256]:
# determining the largest genome
!cd $genomeDirRn; \
    find . -name "*.fna" | \
    seq_tools fasta_info --fn --tn --tl --tgc --sn --sl --sgc --header -n 24 - \
    > genome_info.txt

loading table and plotting


In [260]:
%%R -i genomeDirRn

inFile = file.path(genomeDirRn, 'genome_info.txt')

df = read.delim(inFile, sep='\t')
df %>% head(n=3)


                                   file_name total_sequences total_seq_length
1      ./Streptococcus_pyogenes_MGAS8232.fna               1          1895017
2         ./Helicobacter_cinaedi_PAGU611.fna               1          2078348
3 ./Sphaerobacter_thermophilus_DSM_20745.fna               2          3993764
  total_GC
1    38.55
2    38.63
3    68.11
                                                                       seq_name
1    NC_003485_Streptococcus_pyogenes_MGAS8232__Streptococcus_pyogenes_MGAS8232
2  NC_017761_Helicobacter_cinaedi_PAGU611_DNA__Helicobacter_cinaedi_PAGU611_DNA
3 NC_013523_Sphaerobacter_thermophilus_DSM_20745_chromosome_1__complete_sequenc
  seq_length seq_GC
1    1895017  38.55
2    2078348  38.63
3    2741033  68.10

In [263]:
%%R -w 500 -h 350
# genome lengths
ggplot(df, aes(total_seq_length)) +
    geom_histogram(binwidth=50000) +
    labs(x='Genome length (bp)', y='Count') +
    theme_bw() +
    theme(
        text = element_text(size=16)
    )



In [264]:
%%R
summary(df$total_seq_length)


    Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
  709800  2449000  3732000  3957000  4989000 14780000 

In [265]:
%%R -w 500 -h 350
# genome GC
ggplot(df, aes(total_GC)) +
    geom_histogram(binwidth=1) +
    labs(x='Genome G+C', y='Count') +
    theme_bw() +
    theme(
        text = element_text(size=16)
    )



In [266]:
%%R -w 500 -h 350
# genome GC
ggplot(df, aes(total_GC)) +
    geom_density(fill='red', alpha=0.3) +
    scale_x_continuous(limits=c(10,90)) +
    labs(x='Genome G+C', y='Density') +
    theme_bw() +
    theme(
        text = element_text(size=16)
    )



In [268]:
%%R -w 500 -h 350
# sequences per genome
ggplot(df, aes(total_sequences)) +
    geom_histogram(binwidth=1) +
    labs(x='Total chromosomes/contigs', y='Count') +
    theme_bw() +
    theme(
        text = element_text(size=16)
    )


Indexing genomes


In [271]:
# list of all genomes files and their associated names
!cd $genomeDirRn; \
    find . -name "*fna" | \
    perl -pe 's/.+\///' | \
    perl -pe 's/(.+)(\.[^.]+)/\$1\t\$1\$2/' > genome_index.txt

In [ ]:
!cd $genomeDirRn; \
    SIPSim genome_index \
    genome_index.txt \
    --fp . --np 30 > index_log.txt

Summary


In [ ]:
!printf 'Genome directory: '; echo $genomeDirRn

!printf 'Number of genomes: '
!find $genomeDirRn -name "*.fna"| wc -l

!printf 'Number of indexed genomes: '
!find $genomeDirRn -name "*.fna.sqlite3.db*" | wc -l



-- OLD --

Weighted sampling of genome dataset

  • weighting against low GC organisms

In [285]:
%%R

curve(-exp(-x/10 + 1), from = 10, to = 70)



In [293]:
%%R -h 300
#weights = sapply(df.bac.complete.rep$gc, function(gc) ifelse(gc < 35, 1 - ((35 - gc) / 100), 1))
#weights = sapply(df.bac.complete.rep$gc, function(gc) ifelse(gc < 35, 1 - exp(-gc/10), 1))
#weights = sapply(df.bac.complete.rep$gc, function(gc) 1 - exp(-gc/50 + 2))
weights = data.frame('weights' = weights, 'gc' = df.bac.complete.rep$gc)
    
ggplot(weights, aes(gc, weights)) +
    geom_point() +
    theme_bw()



In [244]:
%%R -w 600 -h 300
n.sample = 500

#weights = sapply(df.bac.complete.rep$gc, function(gc) ifelse(gc < 35, 1 - ((35 - gc) / 50), 1))
to.keep = sample(df.bac.complete.rep$name, n.sample, replace=FALSE, prob=weights$weights)
    
df.bac.complete.rep.p = df.bac.complete.rep %>%
    filter(name %in% to.keep)
    
# plotting GC distribution
ggplot(df.bac.complete.rep.p, aes(gc)) +
    geom_histogram(binwidth=1) +
    geom_vline(xintercept=50, linetype='dashed', color='red', alpha=0.7) +
    labs(x='G+C') +
    theme(
        text = element_text(size=16)
    )



In [ ]: