<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 Description
Simulating CsCl gradient to reproduce the results in
Step1:  Init
Step2:  Downloading genomes
M. barkeri DSM 800
RefSeq = NC_007355.1, NC_007349.1
Methylobacterium extorquens AM1 DSM 1338
RefSeq = NC_012808.1, NC_012811.1, NC_012807.1, NC_012809.1, NC_012810.1
Step3:  Indexing genomes
Step4:  Simulating a gradient community
Step6:  Simulating isotope incorporation
Step7:  Resimulating community (n=2)
Step8:  Simulating gradient fractions
Step9:  Simulating shotgun fragments & calculating GC
Step10:  Creating OTU table
Step11:  Plotting OTU abundances
<ASSISTANT_TASK:>
Python Code:
workDir = "/home/nick/notebook/SIPSim/t/M.bark_M.ext/"
Explanation: Description
Simulating CsCl gradient to reproduce the results in:
Lueders T, Manefield M, Friedrich MW. (2004). Enhanced sensitivity of DNA- and rRNA-based stable isotope probing by fractionation and quantitative analysis of isopycnic centrifugation gradients. Environmental Microbiology 6:73–78.
rotor:
TV865
rotor speed:
45000 rpm
spin time:
>36 hr
gradient average density:
1.725 g/ml
Fraction sizes:
400 ul
Used gDNA
13C-labeled
13C-methanol
5 ug
strains
Methylobacterium extorquens AM1 DSM 1338
M. barkeri DSM 800
Setting variables
End of explanation
import os
import sys
%load_ext rpy2.ipython
%%R
library(ggplot2)
library(dplyr)
library(tidyr)
Explanation: Init
End of explanation
!cd $workDir; \
    seqDB_tools accession-GI2fasta < M.barkeri_refseq.txt > M.barkeri.fna
!cd $workDir; \
    seqDB_tools accession-GI2fasta < M.extorquens_AM1_refseq.txt > M.extorquens_AM1.fna
Explanation: Downloading genomes
M. barkeri DSM 800
RefSeq = NC_007355.1, NC_007349.1
Methylobacterium extorquens AM1 DSM 1338
RefSeq = NC_012808.1, NC_012811.1, NC_012807.1, NC_012809.1, NC_012810.1
End of explanation
# renaming genome sequences
!cd $workDir; \
    find . -name "*_rn.fna" |\
    xargs -I % rm -f %
    
!cd $workDir; \
    find . -name "*.fna" |\
    perl -pe 's/\.fna$//' | \
    xargs -P 2 -I % bash -c \
    "SIPSim renameGenomes %.fna > %_rn.fna"
# list of all genomes files and their associated names
!cd $workDir; \
    find . -name "*_rn.fna" | \
    perl -pe 's/.+\///' | \
    perl -pe 's/(.+)(\.[^.]+)/\$1\t\$1\$2/' > genomes_all_list.txt
    
!cd $workDir; head genomes_all_list.txt
!cd $workDir; \
    SIPSim indexGenomes genomes_all_list.txt \
    --np 2 > index_log.txt
Explanation: Indexing genomes
End of explanation
!cd $workDir; \
    SIPSim gradientComms \
    --n_comm 1 \
    --abund_dist uniform \
    --abund_dist_p low:1,high:1 \
    genomes_all_list.txt > comm-n1-unif.txt
    
!cd $workDir; tail comm-n1-unif.txt
Explanation: Simulating a gradient community
End of explanation
# making config file
config = 
[library 1]
  # baseline: no incorp
  
  [[intraPopDist 1]]
  distribution = uniform
    [[[start]]]
      [[[[interPopDist 1]]]]
        distribution = uniform
        start = 0
        end = 0
    [[[end]]]
      [[[[interPopDist 1]]]]
        distribution = uniform
        start = 0
        end = 0
outfile = os.path.join(workDir, 'incorp.config')
outf = open(outfile, 'wb')
outf.write(config)
outf.close()
!cd $workDir; \
    SIPSim isoIncorp \
    comm-n1-unif.txt incorp.config \
    > incorp-n1-unif.txt
    
!cd $workDir; head incorp-n1-unif.txt
%%bash -s "$workDir"
# adding incorp 100% library for N.extorquens
cd $1
cat incorp-n1-unif.txt \
    <(tail -n +2 incorp-n1-unif.txt | \
    perl -pe 's/^1/2/' | \
    perl -pe 's/0\.0/100.0/ if /M.extorquens/') \
    > incorp-n2-unif.txt
    
    
head incorp-n2-unif.txt
Explanation: Simulating isotope incorporation
End of explanation
!cd $workDir; \
    SIPSim gradientComms \
    --n_comm 2 \
    --abund_dist uniform \
    --abund_dist_p low:1,high:1 \
    genomes_all_list.txt > comm-n2-unif.txt
    
!cd $workDir; tail comm-n2-unif.txt
Explanation: Resimulating community (n=2)
End of explanation
!cd $workDir; \
    SIPSim fractions \
    comm-n2-unif.txt \
    > fracs-n2-unif.txt
    
!cd $workDir; head fracs-n2-unif.txt
Explanation: Simulating gradient fractions
End of explanation
!cd $workDir;\
    SIPSim fragGC \
    genomes_all_list.txt \
    --flr 500,None \
    --fld skewed-normal,9000,2500,-5 \
    --nf 50x \
    --np 24 \
    2> shotFragGC_skewN90-25-n5-nS.log \
    > shotFragGC_skewN90-25-n5-nS.pkl
    
!cd $workDir; head shotFragGC_skewN90-25-n5-nS.log
Explanation: Simulating shotgun fragments & calculating GC
End of explanation
!cd $workDir; \
    SIPSim OTU_sim \
    shotFragGC_skewN90-25-n5-nS.pkl \
    comm-n2-unif.txt \
    incorp-n2-unif.txt \
    fracs-n2-unif.txt \
    --abs_abund 2e8 \
    2> OTU-n2-unif_skewN90-25-n5-nS_A2e9.log \
    > OTU-n2-unif_skewN90-25-n5-nS_A2e9.txt
    
!cd $workDir; head OTU-n2-unif_skewN90-25-n5-nS_A2e9.log 
Explanation: Creating OTU table
End of explanation
%%R -i workDir
# loading file
inFiles = c('OTU-n2-unif_skewN90-25-n5-nS_A2e9.txt')
inFiles = sapply(inFiles, function(x){
            x = as.character(x)
            paste(c(workDir, x), collapse='/')
    })
tbls = list()
for (fileName in inFiles){
    tbls[[fileName]] = read.csv(fileName, sep='\t')
    }
tbl = do.call(rbind, tbls)
tbl$abs_abund = as.numeric(gsub('.+-nS_A|\\.txt\\.[0-9]+', '', rownames(tbl)))
tbl = tbl %>% 
    filter(!grepl('inf', fractions, ignore.case=T)) %>%
    separate(fractions, into = c('BD_min','BD_max'), sep='-', convert=TRUE) %>%
    filter(BD_min != 1.795)
%%R
## BD for G+C of 0 or 100
BD.GCp0 = 0 * 0.098 + 1.66
BD.GCp100 = 1 * 0.098 + 1.66
%%R -w 800 -h 400
# plotting absolute abundances
tbl.s = tbl %>%
    mutate(BD_mean = (BD_min + BD_max) / 2) %>%
    group_by(abs_abund, BD_mean, library, taxon) %>%
    summarize(total_count = sum(count))
## plot
p = ggplot(tbl.s, aes(BD_mean, total_count, shape=taxon, color=taxon)) +
    geom_point() +
    geom_line() +
    scale_x_continuous(limits=c(1.68,1.78), breaks=seq(1.68,1.78,0.02)) +
    labs(x='Buoyant density') +
    facet_grid(library ~ .) +
    theme( 
        text = element_text(size=16) 
    )
p
Explanation: Plotting OTU abundances
End of explanation
<END_TASK>