<jupyter_start><jupyter_text>Description:* Demultiplexing libraries from MiSeq run  Setting variables (user input)<jupyter_code>import os

seqdir = '/home/backup_files/raw_reads/hempmicrobiome.Sam.Ali.SmartLab.2018/'
workdir = '/home/sam/notebooks/hemp_microbiome/data/OTU_pipeline_output/'

# metadata file location(s)
metadata = 'metadata_16S.txt'

index_read1_file = 'index1.16S.fq.gz'
index_read2_file = 'index2.16S.fq.gz'
read12_screed = 'pear_merged-2018-03-30.assembled.fastq_screed'
index_read1_screed = 'index1.16S.fq_screed'
index_read2_screed = 'index1.16S.fq_screed'<jupyter_output><empty_output><jupyter_text>Init<jupyter_code>import os
import json
import screed
import pandas as pd
from collections import Counter
import pandas as pd
import matplotlib.pyplot as plt
import numpy as np
import re
import scipy.stats as stats
import screed
from glob import glob
%cd $workdir
%load_ext rpy2.ipython
#%load_ext pushnote 
%%R
library(ggplot2)
library(dplyr)
library(tidyr)<jupyter_output>/home/sam/anaconda3/lib/python3.6/site-packages/rpy2/rinterface/__init__.py:145: RRuntimeWarning: 
Attaching package: ‘dplyr’


  warnings.warn(x, RRuntimeWarning)
/home/sam/anaconda3/lib/python3.6/site-packages/rpy2/rinterface/__init__.py:145: RRuntimeWarning: The following objects are masked from ‘package:stats’:

    filter, lag


  warnings.warn(x, RRuntimeWarning)
/home/sam/anaconda3/lib/python3.6/site-packages/rpy2/rinterface/__init__.py:145: RRuntimeWarning: The following objects are masked from ‘package:base’:

    intersect, setdiff, setequal, union


  warnings.warn(x, RRuntimeWarning)<jupyter_text>Loading metadata tables<jupyter_code># loading metadata sheet
df_meta = pd.read_csv(workdir+metadata, delimiter="\t")
print(df_meta.head(n=3))<jupyter_output>sample_ID  Extraction_plate Extraction_plate_well  Dilution_plate  \
0    1-2-RT  hemp_DNAx_plate1                    A1  hemp_DI_plate1   
1     6-5-F  hemp_DNAx_plate1                    B1  hemp_DI_plate1   
2     2-3-L  hemp_DNAx_plate1                    C1  hemp_DI_plate1   

  Dilution_plate_well  PCR_run  PCR_order  Primer        Library     Exp_type  \
0                  A1        1          1       1  hemp_library1  whole_plant   
1                  B1        1          2       2  hemp_library1  whole_plant   
2                  C1        1          3       3  hemp_library1  whole_plant   

         Field  Plant  Sample_type  Week  Sex  DNA_conc  Date_extracted  \
0      McGowan    2.0  root_tissue   NaN  NaN     1.712          180111   
1    Freeville    5.0      flowers   NaN  NaN     0.065          180111   
2  East_Ithaca    3.0       leaves   NaN  NaN     0.126          180111   

   DI_factor    Index1    Index2  
0      1.712  ATCGTACG  AACTCTCG  
1      1.000  ACTAT[...]<jupyter_text>Uncompress index files and create screed db<jupyter_code>def uncomp(fileName, directory):
    output = fileName.replace(".gz", "")
    fileName = directory+fileName
    !cd $workdir; \
        pigz -k -d -p 24 -c -f $fileName > $output
    return os.path.splitext(output)[0]
    
index_read1_file = uncomp(index_read1_file, seqdir)
index_read2_file = uncomp(index_read2_file, seqdir)
os.chdir(workdir)
screed.read_fastq_sequences(index_read1_file+".fq")
screed.read_fastq_sequences(index_read2_file+".fq")

index_read1_screed = index_read1_file + '.fq_screed'
index_read2_screed = index_read2_file + '.fq_screed'<jupyter_output><empty_output><jupyter_text>Demultiplex<jupyter_code># rev-comp functions
def complement(seq): 
    idx = {'A': 'T', 'C': 'G', 'G': 'C', 'T': 'A', 'N':'N'} 
    idx.update( {x.lower():y.lower() for x,y in idx.items()} )
    return ''.join(idx[x] for x in list(seq))

def revComp(seq):
    seq = seq[::-1]
    return complement(seq)    
# loading/editing metadata

df_meta.columns = ['Sample'] + list((df_meta.columns)[1:])

df_meta['Barcode_FR'] = [(revComp(y) + x).lower() for x,y \
                         in zip(df_meta.Index1,df_meta.Index2)]


## format of barcode=>sample dict: [revBarcode=read1 + fwdBarcode=read2] : sampleID
map_d = {x:y for x,y in zip(df_meta.Barcode_FR, df_meta.Sample)}
# resulting dict of sample => barcode_FR 
n_print = 96
[(list(map_d.keys())[i], map_d[list(map_d.keys())[i]]) for i in range(n_print)]
# loading screed databases 
ir1db = screed.ScreedDB(index_read1_screed)
ir2db = screed.ScreedDB(index_read2_screed)
amp_db = screed.ScreedDB(read12_screed)
# demultiplexing: naming reads by sample

outFile = re.sub('fastq_screed', 'dmult.fastq', read12_screed)

counter = 0
unassigned = 0
cnt = Counter()
outFile = os.path.join(workdir, os.path.split(outFile)[1])
c = 0
with open(outFile, "w") as out:
    for rec in amp_db.itervalues():
        index_concat = ir1db[rec["name"]]["sequence"] + ir2db[rec["name"]]["sequence"]
        
        try:
            new_name = map_d[index_concat.lower()] + "_" + str(counter) #case sensitive
        except KeyError:
            unassigned += 1
            continue
        counter += 1
        cnt[map_d[index_concat.lower()]] += 1 #case sensitive
        s, q = rec["sequence"], rec["quality"]
        out.write("@%s orig_name=%s\n%s\n+\n%s\n"%(new_name,rec["name"],s,q))
        
        
print ("Unable to assign %s reads to samples"%unassigned)
# number of sequences
n_lines = !wc -l $outFile
n_lines = int(re.sub(' .+', '', n_lines[0]))
print('Number of sequences: {}'.format(n_lines/4))<jupyter_output>Number of sequences: 12638894.0<jupyter_text>Stats on sequences  Number of sequences per sample<jupyter_code># counting sequences for each sample
re_seqName = re.compile(r'_\d+ orig_name.+')

seq_cnt = dict()
with open(outFile, 'r') as fn:
    for line in fn.readlines():
        if re.search(re_seqName, line):
            sampleName = re.sub(re_seqName, '', line).rstrip().lstrip('@')
            try:
                seq_cnt[sampleName] += 1
            except KeyError:
                seq_cnt[sampleName] = 1 
# converting to dataframe
df_seq_cnt = pd.DataFrame.from_dict(seq_cnt, orient='index')
df_seq_cnt.index.name = 'Sample'
df_seq_cnt.reset_index(inplace=True)
df_seq_cnt.columns = ['Sample', 'seq_count']
df_seq_cnt
%%R -i df_seq_cnt -w 900 -h 350

df_seq_cnt$Sample = reorder(df_seq_cnt$Sample, desc(df_seq_cnt$seq_count))

ggplot(df_seq_cnt, aes(Sample, seq_count)) +
    geom_bar(stat='identity') +
    theme_bw() +
    theme( 
        text = element_text(size=16),
        axis.text.x = element_text(angle=45, hjust=1)
    )
%%R
df_seq_cnt[df_seq_cnt$seq_count < 10000,]<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.11968697253337425, 'fraction_numerical': 0.0322234156820623, 'mean_word_length': 3.446111869031378, 'pattern_counts': {'":': 0, '<': 23, '<?xml version=': 0, '>': 25, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 6, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '6', 'dir': 'jupyter-structured-clean-dedup', 'id': '499627', 'license': 'MIT', 'n_tokens_mistral': 2549, 'n_tokens_neox': 2417, 'n_words': 482, 'path': '16S_sequence_processing/2_demultiplex.ipynb', 'repo_name': 'seb369/hemp_microbiome', 'size': '6517'}