<jupyter_start><jupyter_text>Goal: * ID and pull out the sequences of all 16S genes in bac_genome_n1210 dataset  User variables<jupyter_code>import os
baseDir = '/home/nick/notebook/SIPSim/dev/bac_genome1147/'
genomeDir = '/var/seq_data/ncbi_db/genome/Jan2016/bac_complete_spec-rep1_rn/'
rnammerDir = os.path.join(baseDir, 'rnammer')<jupyter_output><empty_output><jupyter_text>Init<jupyter_code>import glob
from IPython.display import Image
%load_ext rpy2.ipython
%%R
library(ggplot2)
library(dplyr)
library(tidyr)
library(gridExtra)
if not os.path.isdir(rnammerDir):
    os.makedirs(rnammerDir)<jupyter_output><empty_output><jupyter_text>rnammer run<jupyter_code>%%bash -s "$genomeDir" "$rnammerDir"

find $1 -name "*.fna" | \
    perl -pe 's/.+\/|\.fna//g' | \
    xargs -n 1 -I % -P 30 bash -c \
    "rnammer -S bac -m ssu -gff $2/%_rrn.gff -f $2/%_rrn.fna -xml $2/%_rrn.xml < $1/%.fna"
## Summarizing the results

!cd $rnammerDir; \
    egrep -v "^#" *.gff | \
    grep "16s_rRNA" | \
    perl -pe 's/:/\t/' > ssu_summary.txt
%%bash -s "$rnammerDir"
cd $1

printf "ssu gene length distribution:\n"
cut -d $'\t' -f 7 ssu_summary.txt | NY_misc_perl stats_descriptive<jupyter_output>ssu gene length distribution:
1	min	589.20
1	Q1	1882.60
1	mean	1903.58
1	median	1940.60
1	Q3	1979.00
1	max	2090.20
1	stdev	138.74<jupyter_text>Compiling 16S sequences<jupyter_code>! cd $rnammerDir; \
    cat *_rrn.fna > bac_genome1147_16S.fna
    
!printf "Number of sequences: "
! cd $rnammerDir; \
    grep -c ">" bac_genome1147_16S.fna<jupyter_output>Number of sequences: 4498
StarcoderdataJupyterStructuredCleanDedup