<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 Pipeling to Process Raw Sequences into Phyloseq Object with DADA2
Prep for Import to QIIME2  (Combine two index files)
Import to QIIME2
Demultiplex
Denoise and Merge
Prepare OTU Tables and Rep Sequences  (Note
Step1:  Step 2
Step2:  Step 3
Step3:  Step 4
Step4:  Step 5
Step5:  Step 6
Step6:  Step 7
Step7:  Step 8
Step8:  Step 9
Step9:  Step 10
Step10:  Step 11
Step11:  Step 13
Step12:  Step 14
Step13:  Step 15
<ASSISTANT_TASK:>
Python Code:
import os, re
# Provide the directory for your index and read files (you can do multiple independently in one go)
bioblitz = '/home/roli/BioBlitz.2017/SV_based/'
# Prepare an object with the name of the library, the name of the directory object (created above), and the metadatafile name
#datasets = [['name',directory1,'metadata1','domain of life'],['name',directory2,'metadata2','domain of life']]
datasets = [['bioblitz',bioblitz,'metadata.tsv','bacteria']]
# Ensure your reads files are named accordingly (or modify to suit your needs)
readFile1 = 'read1.fq.gz'
readFile2 = 'read2.fq.gz'
indexFile1 = 'index_read1.fq.gz'
indexFile2 = 'index_read2.fq.gz'
## Enter Minimum Support for Keeping QIIME Classification
# Note: Classifications that do not meet this criteria will simply be retained, but labeled 'putative'
min_support = 0.8
Explanation: Pipeling to Process Raw Sequences into Phyloseq Object with DADA2
Prep for Import to QIIME2  (Combine two index files)
Import to QIIME2
Demultiplex
Denoise and Merge
Prepare OTU Tables and Rep Sequences  (Note: sample names starting with a digit will break this step)
Classify Seqs
100% Appropriated from the "Atacama Desert Tutorial" for QIIME2
Pipeline can handle both 16S rRNA gene and ITS sequences (in theory)
Tested on 515f and 806r
Tested on ITS1
Commands to Install Dependencies
|| QIIME2 ||
conda create -n qiime2-pipeline --file https://data.qiime2.org/distro/core/qiime2-2017.11-conda-linux-64.txt
source activate qiime2-pipeline
** Note: QIIME2 is still actively in development, and I've noticed frequent new releases. Check for the most up-to-date conda install file https://docs.qiime2.org/2017.11/install/native/#install-qiime-2-within-a-conda-environment
|| Copyrighter rrn Database ||
The script will automatically install the curated GreenGenes rrn attribute database
https://github.com/fangly/AmpliCopyrighter
|| rpy2 (don't use conda version) ||
pip install rpy2  
|| phyloseq ||
conda install -c r r-igraph 
Rscript -e "source('http://bioconductor.org/biocLite.R');biocLite('phyloseq')" 
|| R packages ||
ape   (natively installed with in conda environment)
Citations
Caporaso, J. G., Kuczynski, J., Stombaugh, J., Bittinger, K., Bushman, F. D., Costello, E. K., et al. (2010). QIIME allows analysis of high-throughput community sequencing data. Nature methods, 7(5), 335-336.
McMurdie and Holmes (2013) phyloseq: An R Package for Reproducible Interactive Analysis and Graphics of Microbiome Census Data. PLoS ONE. 8(4):e61217
Paradis E., Claude J. & Strimmer K. 2004. APE: analyses of phylogenetics and evolution in R language. Bioinformatics 20: 289-290.
Angly, F. E., Dennis, P. G., Skarshewski, A., Vanwonterghem, I., Hugenholtz, P., & Tyson, G. W. (2014). CopyRighter: a rapid tool for improving the accuracy of microbial community profiles through lineage-specific gene copy number correction. Microbiome, 2(1), 11.
Last Modified by R. Wilhelm on October 12th, 2017
Step 1: User Input
End of explanation
## Note: QIIME takes a single barcode file. The command 'extract_barcodes.py' concatenates the forward and reverse read barcode and attributes it to a single read.
# See http://qiime.org/tutorials/processing_illumina_data.html
for dataset in datasets:
    directory = dataset[1]
    index1 = directory+indexFile1
    index2 = directory+indexFile2
    
    # Run extract_barcodes to merge the two index files
    !python2 /opt/anaconda2/bin/extract_barcodes.py --input_type barcode_paired_end -f $index1 -r $index2 --bc1_len 8 --bc2_len 8 -o $directory/output
    # QIIME2 import requires a directory containing files names: forward.fastq.gz, reverse.fastq.gz and barcodes.fastq.gz 
    !ln -s $directory$readFile1 $directory/output/forward.fastq.gz
    !ln -s $directory$readFile2 $directory/output/reverse.fastq.gz
    
    # Gzip the barcodes files (apparently necessary)
    !pigz -p 5 $directory/output/barcodes.fastq
    # Removed orphaned reads files (not needed)
    !rm $directory/output/reads?.fastq
Explanation: Step 2: Concatenate Barcodes for QIIME2 Pipeline
End of explanation
for dataset in datasets:
    name = dataset[0]
    directory = dataset[1]
    
    os.system(' '.join([
        "qiime tools import",
        "--type EMPPairedEndSequences",
        "--input-path "+directory+"output/",
        "--output-path "+directory+"output/"+name+".qza"
    ]))
    
    # This more direct command is broken by the fact QIIME uses multiple dashes in their arguments (is my theory)
    #!qiime tools import --type EMPPairedEndSequences --input-path $directory/output --output-path $directory/output/$name.qza
     
Explanation: Step 3: Import into QIIME2
End of explanation
########
## Note: The barcode you supply to QIIME is now a concatenation of your forward and reverse barcode.
# Your 'forward' barcode is actually the reverse complement of your reverse barcode and the 'reverse' is your forward barcode. The file 'primers.complete.csv' provides this information corresponding to the Buckley Lab 'primer number'
# This quirk could be corrected in how different sequencing facilities pre-process the output from the sequencer
##
## SLOW STEP (~ 2 - 4 hrs)
##
for dataset in datasets:
    name = dataset[0]
    directory = dataset[1]
    metadata = dataset[2]
    
    os.system(' '.join([
        "qiime demux emp-paired",
        "--m-barcodes-file "+directory+metadata,
        "--m-barcodes-category BarcodeSequence",
        "--i-seqs "+directory+"output/"+name+".qza",
        "--o-per-sample-sequences "+directory+"output/"+name+".demux"
    ]))
    
    # This more direct command is broken by the fact QIIME uses multiple dashes in their arguments (is my theory)
    #!qiime demux emp-paired --m-barcodes-file $directory/$metadata --m-barcodes-category BarcodeSequence --i-seqs $directory/output/$name.qza --o-per-sample-sequences $directory/output/$name.demux
    
Explanation: Step 4: Demultiplex
End of explanation
## Based on the Graph Produced using the Following Command enter the trim and truncate values. Trim refers to the start of a sequence and truncate the total length (i.e. number of bases to remove from end)
# The example in the Atacam Desert Tutorial trims 13 bp from the start of each read and does not remove any bases from the end of the 150 bp reads:
#  --p-trim-left-f 13 \  
#  --p-trim-left-r 13 \
#  --p-trunc-len-f 150 \
#  --p-trunc-len-r 150
for dataset in datasets:
    name = dataset[0]
    directory = dataset[1]
    
    os.system(' '.join([
        "qiime demux summarize",
        "--i-data "+directory+"/output/"+name+".demux.qza",
        "--o-visualization "+directory+"/output/"+name+".demux.QC.summary.qzv"
    ]))
    
    ## Take the output from this command and drop it into:
    #https://view.qiime2.org
wait_for_user = input("The script will now wait for you to input trimming parameters in the next cell. You will need to take the .qzv files for each library and visualize them at <https://view.qiime2.org>. This is hopefully temporary, while QIIME2 developers improve on q2view.\n\n[ENTER ANYTHING. THIS IS ONLY MEANT TO PAUSE THE PIPELING]")
print("\nThe script is now proceeding. Stay tuned to make sure trimming works.")
Explanation: Step 5: Visualize Quality Scores and Determine Trimming Parameters
End of explanation
## User Input Required
trim_dict = {}
## Input your trimming parameters into a python dictionary for all libraries
#trim_dict["LibraryName1"] = [trim_forward, truncate_forward, trim_reverse, truncate_reverse]
#trim_dict["LibraryName2"] = [trim_forward, truncate_forward, trim_reverse, truncate_reverse]
## Example
trim_dict["bioblitz"] = [1, 240, 1, 190]
Explanation: Step 6: Trimming Parameters | USER INPUT REQUIRED
End of explanation
## Hack for Multithreading
# I hardcoded 'nthreads' in both versions of 'run_dada_paired.R' (find your versions by running 'locate run_dada_paired.R' from your home directory)
# I used ~ 20 threads and the processing finished in ~ 7 - 8hrs
##
## SLOW STEP (~ 6 - 8 hrs, IF multithreading is used)
##
for dataset in datasets:
    name = dataset[0]
    directory = dataset[1]
    
    os.system(' '.join([
        "qiime dada2 denoise-paired",
        "--i-demultiplexed-seqs "+directory+"/output/"+name+".demux.qza",
        "--o-table "+directory+"/output/"+name+".table",
        "--o-representative-sequences "+directory+"/output/"+name+".rep.seqs.final",
        "--p-trim-left-f "+str(trim_dict[name][0]),
        "--p-trim-left-r "+str(trim_dict[name][2]),
        "--p-trunc-len-f "+str(trim_dict[name][1]),
        "--p-trunc-len-r "+str(trim_dict[name][3])
    ]))
  
Explanation: Step 7: Trim, Denoise and Join (aka 'Merge') Reads Using DADA2
End of explanation
for dataset in datasets:
    name = dataset[0]
    directory = dataset[1]
    metadata = dataset[2]
    
    os.system(' '.join([
        "qiime feature-table summarize",
        "--i-table "+directory+"/output/"+name+".table.qza",
        "--o-visualization "+directory+"/output/"+name+".table.qzv",
        "--m-sample-metadata-file "+directory+metadata
    ]))
    os.system(' '.join([
        "qiime feature-table tabulate-seqs",
        "--i-data "+directory+"/output/"+name+".rep.seqs.final.qza",
        "--o-visualization "+directory+"/output/"+name+".rep.seqs.final.qzv"
    ])) 
Explanation: Step 8: Create Summary of OTUs
End of explanation
## Hack for Multithreading
# I hardcoded 'n_threads' in '_mafft.py' in the directory ~/anaconda3/envs/qiime2-2017.9/lib/python3.5/site-packages/q2_alignment
# I used ~ 20 threads and the processing finished in ~ 15 min
for dataset in datasets:
    name = dataset[0]
    directory = dataset[1]
    metadata = dataset[2]
    domain = dataset[3]
    if domain != "fungi":
        # Generate Alignment with MAFFT
        os.system(' '.join([
            "qiime alignment mafft",
            "--i-sequences "+directory+"/output/"+name+".rep.seqs.final.qza",
            "--o-alignment "+directory+"/output/"+name+".rep.seqs.aligned.qza"
        ]))
        # Mask Hypervariable parts of Alignment
        os.system(' '.join([
            "qiime alignment mask",
            "--i-alignment "+directory+"/output/"+name+".rep.seqs.aligned.qza",
            "--o-masked-alignment "+directory+"/output/"+name+".rep.seqs.aligned.masked.qza"
        ])) 
        # Generate Tree with FastTree
        os.system(' '.join([
            "qiime phylogeny fasttree",
            "--i-alignment "+directory+"/output/"+name+".rep.seqs.aligned.masked.qza",
            "--o-tree "+directory+"/output/"+name+".rep.seqs.tree.unrooted.qza"
        ])) 
        # Root Tree
        os.system(' '.join([
            "qiime phylogeny midpoint-root",
            "--i-tree "+directory+"/output/"+name+".rep.seqs.tree.unrooted.qza",
            "--o-rooted-tree "+directory+"/output/"+name+".rep.seqs.tree.final.qza"
        ])) 
Explanation: Step 9: Make Phylogenetic Tree
End of explanation
for dataset in datasets:
    name = dataset[0]
    directory = dataset[1]
    metadata = dataset[2]
    domain = dataset[3]
    # Classify
    if domain == 'bacteria':
        os.system(' '.join([
            "qiime feature-classifier classify-sklearn",
            "--i-classifier /home/db/GreenGenes/qiime2_13.8.99_515.806_nb.classifier.qza",
            "--i-reads "+directory+"/output/"+name+".rep.seqs.final.qza",
            "--o-classification "+directory+"/output/"+name+".taxonomy.final.qza"
        ]))
    if domain == 'fungi':
        os.system(' '.join([
            "qiime feature-classifier classify-sklearn",
            "--i-classifier /home/db/UNITE/qiime2_unite_ver7.99_20.11.2016_classifier.qza",
            "--i-reads "+directory+"/output/"+name+".rep.seqs.final.qza",
            "--o-classification "+directory+"/output/"+name+".taxonomy.final.qza"
        ]))
    # Output Summary
    os.system(' '.join([
        "qiime metadata tabulate",
        "--m-input-file "+directory+"/output/"+name+".taxonomy.final.qza",
        "--o-visualization "+directory+"/output/"+name+".taxonomy.final.summary.qzv"
    ])) 
Explanation: Step 10: Classify Seqs
End of explanation
## Make Function to Re-Format Taxonomy File to Contain Full Column Information 
# and factor in the certain of the taxonomic assignment
def format_taxonomy(tax_file, min_support):
    output = open(re.sub(".tsv",".fixed.tsv",tax_file), "w")
    output.write("\t".join(["OTU","Domain","Phylum","Class","Order","Family","Genus","Species"])+"\n")
    
    with open(tax_file, "r") as f:
        next(f) #skip header
        for line in f:
            line = line.strip()
            line = line.split("\t")
            read_id = line[0]
            tax_string = line[1]
            # Annotate those strings which do not meet minimum support
            if float(line[2]) < float(min_support):
                tax_string = re.sub("__","__putative ",tax_string)
            # Remove All Underscore Garbage (gimmie aesthetics)
            tax_string = re.sub("k__|p__|c__|o__|f__|g__|s__","",tax_string) 
            # Add in columns containing unclassified taxonomic information
            # Predicated on maximum 7 ranks (Domain -> Species)
            full_rank = tax_string.split(";")
            last_classified = full_rank[len(full_rank)-1]
            count = 1
            while last_classified == " ":
                last_classified = full_rank[len(full_rank)-count]
                count = count + 1
            for n in range(full_rank.index(last_classified)+1, 7, 1):
                try:
                    full_rank[n] = "unclassifed "+last_classified
                except:
                    full_rank.append("unclassifed "+last_classified)
            output.write(read_id+"\t"+'\t'.join(full_rank)+"\n")
            
    return()
#####################
## Export from QIIME2
for dataset in datasets:
    name = dataset[0]
    directory = dataset[1]
    metadata = dataset[2]
    domain = dataset[3]
    ## Final Output Names
    fasta_file = directory+"/output/"+name+".rep.seqs.final.fasta"
    tree_file = directory+"/output/"+name+".tree.final.nwk"
    tax_file = directory+"/output/"+name+".taxonomy.final.tsv"
    count_table = directory+"/output/"+name+".counts.final.biom"
    # Export Classifications
    os.system(' '.join([
        "qiime tools export",
        directory+"/output/"+name+".taxonomy.final.qza",
        "--output-dir "+directory+"/output/"
    ]))
    
    # Reformat Classifications to meet phyloseq format
    format_taxonomy(directory+"/output/taxonomy.tsv", min_support)
    # Export SV Table
    os.system(' '.join([
        "qiime tools export",
        directory+"/output/"+name+".table.qza",
        "--output-dir "+directory+"/output/"
    ]))
    # Export SV Sequences
    os.system(' '.join([
        "qiime tools export",
        directory+"/output/"+name+".rep.seqs.final.qza",
        "--output-dir "+directory+"/output/"
    ]))
    
    # Export Tree
    os.system(' '.join([
        "qiime tools export",
        directory+"/output/"+name+".rep.seqs.tree.final.qza",
        "--output-dir "+directory+"/output/"
    ]))
    
    # Rename Exported Files
    %mv $directory/output/dna-sequences.fasta $fasta_file
    %mv $directory/output/feature-table.biom $count_table
    %mv $directory/output/taxonomy.fixed.tsv $tax_file
    
    if domain == "bacteria":
        %mv $directory/output/tree.nwk $tree_file
    
Explanation: Step 11: Prepare Data for Import to Phyloseq
End of explanation
## This step is based on the database contructed for the software 'copyrighter'
## The software itself lacked information about datastructure (and, the import of a biom from QIIME2 failed, likely because there are multiple versions of the biom format)
downloaded = "N"
for dataset in datasets:
    name = dataset[0]
    directory = dataset[1]
    domain = dataset[3]
    if domain == 'bacteria':
        if downloaded == "N":
            ## Download copyrighter database
            !git clone https://github.com/fangly/AmpliCopyrighter $directory/temp/
            ## There are multiple GreenGenes ID numbers for a given taxonomic string.
            ## However, the copyrighter database uses the same average rrn copy number.
            ## We will therefore just use the taxonomic strings, since QIIME2 does not output the ID numbers
            !sed -e '1,1075178d; 1078115d' $directory/temp/data/201210/ssu_img40_gg201210.txt > $directory/output/copyrighter.tax.strings.tsv
            ## Create Dictionary of rrnDB
            rrnDB = {}
            with open(directory+"/output/copyrighter.tax.strings.tsv", "r") as f:
                for line in f:
                    line = line.strip()
                    line = line.split("\t")
                    try:
                        rrnDB[line[0]] = line[1]
                    except:
                        pass
            downloaded = "Y"
            
        ## Attribute rrn to readID from taxonomy.tsv
        output = open(directory+"/output/"+name+".seqID.to.rrn.final.tsv","w")
        output.write("Feature ID\trrn\n")
        with open(directory+"/output/taxonomy.tsv", "r") as f:
            missing = 0
            total = 0
            next(f)  # Skip Header
            for line in f:
                line = line.strip()
                line = line.split("\t")
                seqID = line[0]
                try:
                    rrn = rrnDB[line[1]]
                except:
                    rrn = "NA"
                    missing = missing + 1
                total = total + 1
                output.write(seqID+"\t"+rrn+"\n")
        print("\nPercent of OTUs Missing {:.1%}".format(float(missing)/total))
        print("Don't Panic! The majority of missing OTUs could be low abundance.")
Explanation: Step 13: Get 16S rRNA Gene Copy Number (rrn)
End of explanation
## Setup R-Magic for Jupyter Notebooks
import rpy2
%load_ext rpy2.ipython
def fix_biom_conversion(file):
    with open(file, 'r') as fin:
        data = fin.read().splitlines(True)
    with open(file, 'w') as fout:
        fout.writelines(data[1:])
import pandas as pd
%R library(phyloseq)
%R library(ape)
for dataset in datasets:
    name = dataset[0]
    directory = dataset[1]
    metadata = dataset[2]
    domain = dataset[3]
 
    #### IMPORT DATA to R
    ## For '.tsv' files, use Pandas to create a dataframe and then pipe that to R
    ## For '.biom' files, first convert using 'biom convert' on the command-line
    ## Had problems importing the count table with pandas, opted for using read.table in R
    
    # Import Taxonomy File
    tax_file = pd.read_csv(directory+"/output/"+name+".taxonomy.final.tsv", sep="\t")
    %R -i tax_file
    %R rownames(tax_file) = tax_file$OTU
    %R tax_file$OTU <- NULL
    %R tax_file <- tax_file[sort(row.names(tax_file)),] #read names must match the count_table
    # Import Sample Data
    #sample_file = pd.read_csv(directory+"/"+metadata, sep="\t")
    sample_file = pd.read_table(directory+metadata, keep_default_na=False)
    %R -i sample_file
    %R rownames(sample_file) = sample_file$X.SampleID   
    %R sample_file$X.SampleID <- NULL
    %R sample_file$LinkerPrimerSequence <- NULL  ## Clean-up some other stuff
    
    # Import Count Data
    os.system(' '.join([
        "biom convert",
        "-i",
        directory+"/output/"+name+".counts.final.biom",
        "-o",
        directory+"/output/"+name+".counts.final.tsv",
        "--to-tsv"
    ]))
    
    # The biom converter adds a stupid line that messes with the table formatting
    fix_biom_conversion(directory+"/output/"+name+".counts.final.tsv")
    # Finally import
    count_table = pd.read_csv(directory+"/output/"+name+".counts.final.tsv", sep="\t")
    %R -i count_table
    %R rownames(count_table) = count_table$X.OTU.ID   
    %R count_table$X.OTU.ID <- NULL    
    %R count_table <- count_table[sort(row.names(count_table)),] #read names must match the tax_table
    
    # Convert to Phyloseq Objects
    %R p_counts = otu_table(count_table, taxa_are_rows = TRUE)    
    %R p_samples = sample_data(sample_file)    
    %R p_tax = tax_table(tax_file)
    %R taxa_names(p_tax) <- rownames(tax_file) # phyloseq throws out rownames
    %R colnames(p_tax) <- colnames(tax_file) # phyloseq throws out colnames
    
    # Merge Phyloseq Objects
    %R p = phyloseq(p_counts, p_tax)
    # Import Phylogenetic Tree
    if domain == "bacteria":
        tree_file = directory+"/output/"+name+".tree.final.nwk"
        %R -i tree_file  
        %R p_tree <- read.tree(tree_file)
    
        # Combine All Objects into One Phyloseq
        %R p_final <- merge_phyloseq(p, p_samples, p_tree)
    
    else:
        # Combine All Objects into One Phyloseq
        %R p_final <- merge_phyloseq(p, p_samples)
        
    # Save Phyloseq Object as '.rds'
    output = directory+"/output/p_"+name+".final.rds"
    %R -i output
    %R saveRDS(p_final, file = output)
    
    # Confirm Output
    %R print(p_final)
Explanation: Step 14: Import into Phyloseq
End of explanation
for dataset in datasets:
    directory = dataset[1]
    metadata = dataset[2]
    
    # Remove Files
    if domain == "bacteria":
        %rm -r $directory/output/*tree.unrooted.qza 
        %rm -r $directory/output/*aligned.masked.qza 
        
    %rm $directory/output/*.biom 
    %rm -r $directory/temp/
    %rm $directory/output/*barcodes.fastq.gz 
    %rm $directory/output/taxonomy.tsv
    %rm $directory/output/forward.fastq.gz # Just the symlink
    %rm $directory/output/reverse.fastq.gz # Just the symlink
    %rm $directory/output/copyrighter.tax.strings.tsv
    
    # Separate Final Files
    %mkdir $directory/final/    
    %mv $directory/output/*.final.rds $directory/final/
    %mv $directory/output/*.taxonomy.final.tsv $directory/final/    
    %mv $directory/output/*.counts.final.tsv $directory/final/
    %mv $directory/output/*.final.fasta $directory/final/
    %cp $directory$metadata $directory/final/
    %mv $directory/output/*.seqID.to.rrn.final.tsv $directory/final/ 
    %mv $directory/output/*.nwk $directory/final/ 
    
    # Gzip and Move Intermediate Files
    !pigz -p 10 $directory/output/*.qza
    !pigz -p 10 $directory/output/*.qzv
    
    %mv $directory/output/ $directory/intermediate_files
print("Your sequences have been successfully saved to 'final' and 'intermediate_files'")
Explanation: Step 15: Clean-up Intermediate Files and Final Outputs
End of explanation
<END_TASK>