<jupyter_start><jupyter_text>Description* Merging raw reads with PEAR  Setting variables<jupyter_code>## Where are your raw sequences found?
seqdir = '/home/backup_files/raw_reads/hempmicrobiome.Sam.Ali.SmartLab.2018/'

## What directory do you want to work in and keep all subsequent files in?
workdir = '/home/sam/notebooks/hemp_microbiome/data/ITS_OTUs/merge_demult/'

## What is the name of your forward and reverse reads?
readFile1 = 'read1.ITS.fq.gz'
readFile2 = 'read2.ITS.fq.gz'

## What name do you want included for all subsequent files?
name = 'hemp_ITS'<jupyter_output><empty_output><jupyter_text>Init<jupyter_code>import screed
from glob import glob
import matplotlib.pyplot as plt
import numpy as np
from mpld3 import enable_notebook
import screed
import pandas as pd
import os
%matplotlib inline  
%load_ext rpy2.ipython
%%R
library(ggplot2)
library(dplyr)
# Move into working directory and if it doesn't exist, make it
if not os.path.isdir(workdir):
    os.makedirs(workdir)
%cd $workdir<jupyter_output>/home/sam/notebooks/hemp_microbiome/data/ITS_OTUs/merge_demult<jupyter_text>Uncompress the fastq files<jupyter_code>output1 = os.path.join(workdir, "forward.fastq")
seqfile1 = os.path.join(seqdir, readFile1)

!cd $workdir; \
    pigz -k -d -p 20 -c $seqfile1 > $output1
output2 = os.path.join(workdir, "reverse.fastq")
seqfile2 = os.path.join(seqdir, readFile2)

!cd $workdir; \
    pigz -k -d -p 20 -c $seqfile2 > $output2<jupyter_output><empty_output><jupyter_text>Merging<jupyter_code>!cd $workdir; \
    pear -m 600 -j 20 \
    -f forward.fastq \
    -r reverse.fastq \
    -o pear_merged_$name<jupyter_output>____  _____    _    ____ 
|  _ \| ____|  / \  |  _ \
| |_) |  _|   / _ \ | |_) |
|  __/| |___ / ___ \|  _ <
|_|   |_____/_/   \_\_| \_\

PEAR v0.9.10 [May 30, 2016]

Citation - PEAR: a fast and accurate Illumina Paired-End reAd mergeR
Zhang et al (2014) Bioinformatics 30(5): 614-620 | doi:10.1093/bioinformatics/btt593

Forward reads file.................: forward.fastq
Reverse reads file.................: reverse.fastq
PHRED..............................: 33
Using empirical frequencies........: YES
Statistical method.................: OES
Maximum assembly length............: 600
Minimum assembly length............: 50
p-value............................: 0.010000
Quality score threshold (trimming).: 0
Minimum read size after trimming...: 1
Maximal ratio of uncalled bases....: 1.000000
Minimum overlap....................: 10
Scoring method.....................: Scaled score
Threads............................: 20

Allocating memory..................: 200,000,000 bytes
Computing empiric[...]<jupyter_text>Making a screed db of merged reads<jupyter_code>pear_merged_file = !echo "pear_merged_"$name".assembled.fastq"
pear_merged_file = pear_merged_file[0]

os.chdir(workdir)
screed.read_fastq_sequences(pear_merged_file)
pear_merged_file += '_screed'
fqdb = screed.ScreedDB(pear_merged_file)
pear_merged_file
lengths = []
for read in fqdb.itervalues():
    lengths.append((len(read["sequence"])))
fig = plt.figure()
ax = fig.add_subplot(111)
h = ax.hist(np.array(lengths), bins=50)
xl = ax.set_xlabel("Sequence Length, nt")
yl = ax.set_ylabel("Count")
fig.set_size_inches((10,6))
print ('Number of reads: {}'.format(len(lengths)))<jupyter_output>Number of reads: 14495332<jupyter_text>Quality stats on merged reads<jupyter_code>def qualStats(sourceDir, fileName):
    outFile = fileName + '_qualStats'
    !cd $sourceDir; \
        fastx_quality_stats -i $fileName -o $outFile -Q 33
    return outFile
    
qualStatsRes = qualStats(workdir, 'pear_merged_'+name+'.assembled.fastq')
%%R -i workdir -i qualStatsRes

setwd(workdir)

# reading in qual-stats files    
tbl.r12 = read.delim(qualStatsRes, sep='\t')
rownames(tbl.r12) = 1:nrow(tbl.r12)
%%R -w 800 -h 300
# smooth curve on median qual values
ggplot(tbl.r12, aes(x=column, y=med, ymin=Q1, ymax=Q3)) +
    geom_smooth(se=FALSE, method='auto') +
    geom_linerange(alpha=0.3) +
    labs(x='position', y='median quality score') +
    theme_bw() +
    theme( 
        text = element_text(size=16) 
        )<jupyter_output>/home/sam/anaconda3/lib/python3.6/site-packages/rpy2/rinterface/__init__.py:145: RRuntimeWarning: `geom_smooth()` using method = 'loess'

  warnings.warn(x, RRuntimeWarning)
{'fraction_non_alphanumeric': 0.17144188213370604, 'fraction_numerical': 0.03377591427905893, 'mean_word_length': 3.5057712486883528, 'pattern_counts': {'":': 0, '<': 22, '<?xml version=': 0, '>': 23, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '6', 'dir': 'jupyter-structured-clean-dedup', 'id': '581458', 'license': 'MIT', 'n_tokens_mistral': 1668, 'n_tokens_neox': 1529, 'n_words': 353, 'path': 'ITS_sequence_processing/1_read_merge_PEAR.ipynb', 'repo_name': 'seb369/hemp_microbiome', 'size': '4293'}