Description:

  • Plotting histograms of fragment lengths for amplicon- and shotgun-fragments

In [3]:
%load_ext rmagic
%load_ext rpy2.ipython


The rmagic extension is already loaded. To reload it, use:
  %reload_ext rmagic

In [4]:
%%R
require(ggplot2)
require(reshape)
require(data.table)
require(gridExtra)

Loading bacteria data

Loading amplicon data


In [5]:
%%R
# loading amplicon data 
base.dir = "~/notebook/deltaGC_notes/data/amplicon/140411_allBac/amp454-L215-N1k";
files = c(
    'skewnorm_m15k_r4k-20k',
    'skewnorm_m10k_r4k-11.5k'
    )

tbl.l = list()
for(i in files){
    in.file = paste(c(i, "txt"), collapse=".")
    in.file = paste(c(base.dir, in.file), collapse="/")
    tbl.l[[i]] = fread(in.file, header=T, sep="\t")
    tbl.l[[i]]$file = rep(i, nrow(tbl.l[[i]]))
    }


## merging lists
tbl.amp = do.call(rbind, tbl.l)
tbl.amp$data = rep('amplicon\nfragments', nrow(tbl.amp))
tbl.l = list()


Read 64.2% of 1091000 rows
Read 98.1% of 1091000 rows
Read 1091000 rows and 11 (of 11) columns from 0.151 GB file in 00:00:04

Loading shotgun data


In [6]:
%%R
# loading shotgun read data
base.dir = "~/notebook/deltaGC_notes/data/mg-reads/140411_allBac/mg454-L215-N10k";
files = list(
    'skewnorm_m15k_r4k-20k',
    'skewnorm_m10k_r4k-11.5k'
    )

tbl.l = list()
for(i in files){
    in.file = paste(c(i, "txt"), collapse=".")
    in.file = paste(c(base.dir, in.file), collapse="/")
    tbl.l[[i]] = fread(in.file, header=T, sep="\t")
    tbl.l[[i]]$file = rep(i, nrow(tbl.l[[i]]))
    }


## merging lists
tbl.mg = do.call(rbind, tbl.l)
tbl.mg$data = rep('shotgun\nfragments', nrow(tbl.mg))
tbl.l = list()


Read 0.0% of 10910000 rows
Read 6.4% of 10910000 rows
Read 9.1% of 10910000 rows
Read 10.9% of 10910000 rows
Read 12.6% of 10910000 rows
Read 14.6% of 10910000 rows
Read 16.5% of 10910000 rows
Read 18.3% of 10910000 rows
Read 20.2% of 10910000 rows
Read 22.4% of 10910000 rows
Read 24.7% of 10910000 rows
Read 27.0% of 10910000 rows
Read 29.6% of 10910000 rows
Read 32.3% of 10910000 rows
Read 35.0% of 10910000 rows
Read 37.9% of 10910000 rows
Read 41.0% of 10910000 rows
Read 44.2% of 10910000 rows
Read 47.6% of 10910000 rows
Read 51.1% of 10910000 rows
Read 54.9% of 10910000 rows
Read 58.8% of 10910000 rows
Read 62.9% of 10910000 rows
Read 67.2% of 10910000 rows
Read 71.7% of 10910000 rows
Read 76.4% of 10910000 rows
Read 80.7% of 10910000 rows
Read 81.5% of 10910000 rows
Read 86.7% of 10910000 rows
Read 92.2% of 10910000 rows
Read 98.0% of 10910000 rows
Read 10910000 rows and 11 (of 11) columns from 1.523 GB file in 00:00:54

Read 0.0% of 10910000 rows
Read 7.6% of 10910000 rows
Read 15.5% of 10910000 rows
Read 23.4% of 10910000 rows
Read 31.3% of 10910000 rows
Read 39.0% of 10910000 rows
Read 46.7% of 10910000 rows
Read 54.2% of 10910000 rows
Read 61.6% of 10910000 rows
Read 69.0% of 10910000 rows
Read 76.4% of 10910000 rows
Read 84.4% of 10910000 rows
Read 92.6% of 10910000 rows
Read 10910000 rows and 11 (of 11) columns from 1.512 GB file in 00:00:18

In [7]:
%%R
## merging amp & mg tables
tbl.bac = rbind(tbl.amp, tbl.mg)
tbl.amp = vector()  # clearing memory
tbl.mg = vector()

In [8]:
%%R
# editing $file names
tbl.bac[, file:=gsub("skewnorm_m15k_r4k-20k", "Fragment size\ndistribution: 12.5 kb", file)]
tbl.bac[, file:=gsub("skewnorm_m10k_r4k-11.5k", "Fragment size\ndistribution: 6 kb", file)]


                                 genome scaffold           read  read_GC
       1:  Agrobacterium_fabrum_str_C58 AE007869  AE007869__499 52.55814
       2:  Agrobacterium_fabrum_str_C58 AE007869  AE007869__498 52.55814
       3:  Agrobacterium_fabrum_str_C58 AE007869  AE007869__495 52.55814
       4:  Agrobacterium_fabrum_str_C58 AE007869  AE007869__497 52.55814
       5:  Agrobacterium_fabrum_str_C58 AE007869  AE007869__496 52.55814
      ---                                                               
24001996: Salmonella_bongori_NCTC_12419 FR877557 FR877557__2790 57.20930
24001997: Salmonella_bongori_NCTC_12419 FR877557 FR877557__2797 29.76744
24001998: Salmonella_bongori_NCTC_12419 FR877557 FR877557__2796 44.18605
24001999: Salmonella_bongori_NCTC_12419 FR877557 FR877557__2795 50.23256
24002000: Salmonella_bongori_NCTC_12419 FR877557 FR877557__2794 53.95349
          read_buoyant_density read_start read_length fragment_GC
       1:             1.711507     114437         215    53.02019
       2:             1.711507     114437         215    55.41002
       3:             1.711507    5035543         215    52.92066
       4:             1.711507    5035543         215    53.45668
       5:             1.711507     114437         215    56.14006
      ---                                                        
24001996:             1.716065    1863392         215    53.96161
24001997:             1.689172    8744011         215    44.52086
24001998:             1.703302    8668917         215    46.68961
24001999:             1.709228    6849897         215    54.76029
24002000:             1.712874    3975166         215    45.70962
          fragment_buoyant_density fragment_start fragment_length
       1:                 1.711960         114273           12284
       2:                 1.714302         107900           14048
       3:                 1.711862        5024735           14911
       4:                 1.712388        5030946           10053
       5:                 1.715017         107683           12109
      ---                                                        
24001996:                 1.712882        1861101            5730
24001997:                 1.703630        8742874            6543
24001998:                 1.705756        8667460            6102
24001999:                 1.713665        6849320            5882
24002000:                 1.704795        3972279            6130
                                          file                data
       1: Fragment size\ndistribution: 12.5 kb amplicon\nfragments
       2: Fragment size\ndistribution: 12.5 kb amplicon\nfragments
       3: Fragment size\ndistribution: 12.5 kb amplicon\nfragments
       4: Fragment size\ndistribution: 12.5 kb amplicon\nfragments
       5: Fragment size\ndistribution: 12.5 kb amplicon\nfragments
      ---                                                         
24001996:    Fragment size\ndistribution: 6 kb  shotgun\nfragments
24001997:    Fragment size\ndistribution: 6 kb  shotgun\nfragments
24001998:    Fragment size\ndistribution: 6 kb  shotgun\nfragments
24001999:    Fragment size\ndistribution: 6 kb  shotgun\nfragments
24002000:    Fragment size\ndistribution: 6 kb  shotgun\nfragments

Loading archaea data

Loading amplicon data


In [9]:
%%R
base.dir = "~/notebook/deltaGC_notes/data/amplicon/140416_allArc/amp454-L215-N1k";
files = c(
    'skewnorm_m15k_r4k-20k',
    'skewnorm_m10k_r4k-11.5k'
    )

tbl.l = list()
for(i in files){
    in.file = paste(c(i, "txt"), collapse=".")
    in.file = paste(c(base.dir, in.file), collapse="/")
    tbl.l[[i]] = fread(in.file, header=T, sep="\t")
    tbl.l[[i]]$file = rep(i, nrow(tbl.l[[i]]))
    }

## merging lists
tbl.amp = do.call(rbind, tbl.l)
tbl.amp$data = rep('amplicon\nfragments', nrow(tbl.amp))
tbl.l = list()

Loading shotgun read data


In [10]:
%%R
base.dir = "~/notebook/deltaGC_notes/data/mg-reads/140417_allArc/mg454-L215-N10k";
files = list(
    'skewnorm_m15k_r4k-20k',
    'skewnorm_m10k_r4k-11.5k'
    )

tbl.l = list()
for(i in files){
    in.file = paste(c(i, "txt"), collapse=".")
    in.file = paste(c(base.dir, in.file), collapse="/")
    tbl.l[[i]] = fread(in.file, header=T, sep="\t")
    tbl.l[[i]]$file = rep(i, nrow(tbl.l[[i]]))
    }


## merging lists
tbl.mg = do.call(rbind, tbl.l)
tbl.mg$data = rep('shotgun\nfragments', nrow(tbl.mg))
tbl.l = list()

In [11]:
%%R
## merging amp & mg tables
tbl.arc = rbind(tbl.amp, tbl.mg)
tbl.amp = vector()
tbl.mg = vector()

In [12]:
%%R
# editing $file names
tbl.arc[, file:=gsub("skewnorm_m15k_r4k-20k", "Fragment size\ndistribution: 12.5 kb", file)]
tbl.arc[, file:=gsub("skewnorm_m10k_r4k-11.5k", "Fragment size\ndistribution: 6 kb", file)]


                                    genome scaffold           read  read_GC
      1: Haloarcula_marismortui_ATCC_43049 AY596297   AY596297__33 54.41860
      2: Haloarcula_marismortui_ATCC_43049 AY596297  AY596297__791 63.72093
      3: Haloarcula_marismortui_ATCC_43049 AY596297   AY596297__37 63.72093
      4: Haloarcula_marismortui_ATCC_43049 AY596297   AY596297__36 63.72093
      5: Haloarcula_marismortui_ATCC_43049 AY596297  AY596297__797 54.41860
     ---                                                                   
2221996:         Haloquadratum_walsbyi_C23 FR746099 FR746099__8237 44.65116
2221997:         Haloquadratum_walsbyi_C23 FR746099 FR746099__8230 49.76744
2221998:         Haloquadratum_walsbyi_C23 FR746099 FR746099__8231 43.25581
2221999:         Haloquadratum_walsbyi_C23 FR746099 FR746099__8232 42.32558
2222000:         Haloquadratum_walsbyi_C23 FR746099 FR746099__8233 40.46512
         read_buoyant_density read_start read_length fragment_GC
      1:             1.713330    2215176         215    57.38426
      2:             1.722447    5265667         215    57.17372
      3:             1.722447    5265667         215    57.25134
      4:             1.722447    5265667         215    56.66473
      5:             1.713330    2215176         215    55.90927
     ---                                                        
2221996:             1.703758    1283771         215    45.85610
2221997:             1.708772    1506868         215    52.50473
2221998:             1.702391    3992795         215    45.10173
2221999:             1.701479     323240         215    47.09707
2222000:             1.699656    5308016         215    46.21100
         fragment_buoyant_density fragment_start fragment_length
      1:                 1.716237        2203273           12960
      2:                 1.716030        5261552           11110
      3:                 1.716106        5261559           10798
      4:                 1.715531        5264670           13774
      5:                 1.714791        2211480           12565
     ---                                                        
2221996:                 1.704939        1277976            6588
2221997:                 1.711455        1506629            5809
2221998:                 1.704200        3987706            5308
2221999:                 1.706155         318123            6459
2222000:                 1.705287        5305459            5384
                                         file                data
      1: Fragment size\ndistribution: 12.5 kb amplicon\nfragments
      2: Fragment size\ndistribution: 12.5 kb amplicon\nfragments
      3: Fragment size\ndistribution: 12.5 kb amplicon\nfragments
      4: Fragment size\ndistribution: 12.5 kb amplicon\nfragments
      5: Fragment size\ndistribution: 12.5 kb amplicon\nfragments
     ---                                                         
2221996:    Fragment size\ndistribution: 6 kb  shotgun\nfragments
2221997:    Fragment size\ndistribution: 6 kb  shotgun\nfragments
2221998:    Fragment size\ndistribution: 6 kb  shotgun\nfragments
2221999:    Fragment size\ndistribution: 6 kb  shotgun\nfragments
2222000:    Fragment size\ndistribution: 6 kb  shotgun\nfragments

Merging bacterial & archaeal tables


In [13]:
%%R
# merging tables
tbl.bac[, domain:= 'Bacteria']
tbl.arc[, domain:= 'Archaea']
#tbl = rbind(tbl.bac, tbl.arc)
#tbl$domain = factor(tbl$domain, levels=c('Bacteria', 'Archaea'))
xxx=1

Plotting fragment length distributions


In [15]:
%%R -h 5 -w 9 -u in
# summarizing fragment length distribution

# plot of fragment lengths
make_frag_len_hist = function(x){
    p = ggplot(x, aes(fragment_length)) +
        geom_histogram(binwidth=100) +
        labs(x='Fragment length (bp)', y='Number of fragments') +
        facet_grid(data ~ file, scales='free') +
        theme(
            text = element_text(size=16),
            axis.text.x = element_text(angle=45, color='black'),
            axis.text.y = element_text(color='black')
            )
    return(p)
    }

# plotting bacterial fragments
p.bac = make_frag_len_hist(tbl.bac)
ggsave(p.bac, file="/home/nick/notebook/deltaGC_notes/data/amp-mg/figures/frag-len-dist_bac.tiff", 
       dpi=300, height=5, width=9, units='in')
print(p.bac)



In [16]:
%%R -h 5 -w 9 -u in
# plotting archaeal fragments
p.arc = make_frag_len_hist(tbl.arc)
ggsave(p.arc, file="/home/nick/notebook/deltaGC_notes/data/amp-mg/figures/frag-len-dist_arc.tiff", 
        dpi=300, height=5, width=9, units='in')
print(p.arc)



In [ ]: