<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 Merge Sequencing Libraries from Different Sequencing Runs that have been denoised by DADA2  and Finish with the Standard Processing Pipeline
Merge multiple libraries for the SAME sequencing region (i.e. ITS+ITS... or SSU+SSU...)
Complete remainder of processing pipieline
Step1:  Step 2
Step2:  Step 3
Step3:  Step 4
Step4:  Step 5
Step5:  Step 6
Step6:  Step 7
Step7:  Step 8
<ASSISTANT_TASK:>
Python Code:
import os, re
# Provide the directory where files are located 
directory = '/home/roli/FORESTs_BHAVYA/Combined_Libraries/ITS/'
#directory = '/home/roli/FORESTs_BHAVYA/Combined_Libraries/16S/'
# Provide a list of all the FeatureTables you will merge
# Produced by QIIME2 in STEP 7 (i.e. DADA2 Denoising/Merging/FeatureTable)
qva_files = ['ITS.table.Honnedaga.qza','ITS.table.Woods.qza']
#qva_files = ['SSU.table.Honnedaga.qza','SSU.table.Woods.qza']
# Provide a list of all the Representative Sequences you will merge
# Also produced from STEP 7 (i.e. DADA2 Denoising/Merging/FeatureTable)
seq_files = ['ITS.rep.seqs.Honnedaga.qza','ITS.rep.seqs.Woods.qza']
#seq_files = ['SSU.rep.seqs.Honnedaga.qza','SSU.rep.seqs.Woods.qza']
# Provide a concatenated metadatafile with matching column order
metadata = 'ITS.metadata.combined.tsv'
## 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
domain = 'fungi'   # 'bacteria' | 'fungi'
Explanation: Pipeling to Merge Sequencing Libraries from Different Sequencing Runs that have been denoised by DADA2  and Finish with the Standard Processing Pipeline
Merge multiple libraries for the SAME sequencing region (i.e. ITS+ITS... or SSU+SSU...)
Complete remainder of processing pipieline: prepare Rep Seqs, Classify Seqs, Make Tree, Make Phyloseq Object
Informed by the "Fecal Microbiota Transplant Tutorial" for QIIME2
Commands to Install Dependencies
|| QIIME2 ||
conda create -n qiime2-2017.9 --file https://data.qiime2.org/distro/core/qiime2-2017.9-conda-linux-64.txt  
source activate qiime2-2017.9
|| rpy2 (don't use conda version) ||
pip install rpy2  
|| phyloseq ||
conda install -c r-igraph 
Rscript -e "source('http://bioconductor.org/biocLite.R');biocLite('phyloseq')" 
|| R packages ||
ape   (natively installed with in conda environment)
Last Modified by R. Wilhelm on October 12th, 2017
Step 1: User Input
End of explanation
%mkdir $directory/output
for n in range(0, len(qva_files), 1):
    if n == 0:
        os.system(' '.join([
            "qiime feature-table merge",
            "--i-table1 "+directory+"/"+qva_files[n],
            "--i-table2 "+directory+"/"+qva_files[n+1],
            "--o-merged-table "+directory+"/output/merged.table.final.qza"
        ]))
        os.system(' '.join([
            "qiime feature-table merge-seq-data",
            "--i-data1 "+directory+"/"+seq_files[n],
            "--i-data2 "+directory+"/"+seq_files[n+1],
            "--o-merged-data "+directory+"/output/merged.rep.seqs.final.qza"
        ]))
    
    if n > 0 and n + 1 < len(qva_files):
        os.system(' '.join([
            "qiime feature-table merge",
            "--i-table1 "+directory+"/output/merged.table.final.qza",
            "--i-table2 "+directory+"/"+qva_files[n+1],
            "--o-merged-table "+directory+"/output/merged.table.final.qza"
        ]))
    
        os.system(' '.join([
            "qiime feature-table merge-seq-data",
            "--i-data1 "+directory+"/output/merged.rep.seqs.final.qza",
            "--i-data2 "+directory+"/"+seq_files[n+1],
            "--o-merged-data "+directory+"/output/merged.rep.seqs.final.qza"
        ]))
Explanation: Step 2: Merge Feature Tables using QIIME
End of explanation
os.system(' '.join([
    "qiime feature-table summarize",
    "--i-table "+directory+"/output/merged.table.final.qza",
    "--o-visualization "+directory+"/output/merged.table.final.qzv",
    "--m-sample-metadata-file "+directory+metadata
]))
os.system(' '.join([
    "qiime feature-table tabulate-seqs",
    "--i-data "+directory+"/output/merged.rep.seqs.final.qza",
    "--o-visualization "+directory+"/output/merged.rep.seqs.final.qzv"
])) 
Explanation: Step 3: Create Summary of OTUs
End of explanation
if domain != "fungi":
    # Generate Alignment with MAFFT
    os.system(' '.join([
        "qiime alignment mafft",
        "--i-sequences "+directory+"/output/merged.rep.seqs.final.qza",
        "--o-alignment "+directory+"/output/merged.rep.seqs.aligned.qza"
    ]))
    # Mask Hypervariable parts of Alignment
    os.system(' '.join([
        "qiime alignment mask",
        "--i-alignment "+directory+"/output/merged.rep.seqs.aligned.qza",
        "--o-masked-alignment "+directory+"/output/merged.rep.seqs.aligned.masked.qza"
    ])) 
    # Generate Tree with FastTree
    os.system(' '.join([
        "qiime phylogeny fasttree",
        "--i-alignment "+directory+"/output/merged.rep.seqs.aligned.masked.qza",
        "--o-tree "+directory+"/output/merged.rep.seqs.tree.unrooted.qza"
    ])) 
    # Root Tree
    os.system(' '.join([
        "qiime phylogeny midpoint-root",
        "--i-tree "+directory+"/output/merged.rep.seqs.tree.unrooted.qza",
        "--o-rooted-tree "+directory+"/output/merged.rep.seqs.tree.final.qza"
    ])) 
Explanation: Step 4: Make Phylogenetic Tree
End of explanation
# 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/merged.rep.seqs.final.qza",
        "--o-classification "+directory+"/output/merged.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/merged.rep.seqs.final.qza",
        "--o-classification "+directory+"/output/merged.taxonomy.final.qza"
    ]))
# Output Summary
os.system(' '.join([
    "qiime metadata tabulate",
    "--m-input-file "+directory+"/output/merged.taxonomy.final.qza",
    "--o-visualization "+directory+"/output/merged.taxonomy.final.summary.qzv"
])) 
Explanation: Step 5: 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):
    rank_dict = {'k__':"Domain",'k__':"Domain",}
    output = open(re.sub(".tsv",".fixed.tsv",tax_file), "w")
    output.write("\t".join(["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)
            
            # Manage cases were taxonomic string is empty
            tax_string = re.sub("k__;|p__;|c__;|o__;|f__;|g__;","unclassified;",tax_string) 
            tax_string = re.sub("s__","unclassified",tax_string)  #note: this makes all species unclassified.
        
            # 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]
            
            for n in range(len(full_rank), 7, 1):
                full_rank.append("unclassifed "+last_classified)
                
            output.write(read_id+"\t"+'\t'.join(full_rank)+"\n")
            
    return()
#####################
## Export from QIIME2
## Final Output Names
fasta_file = directory+"/output/merged.rep.seqs.final.fasta"
tree_file = directory+"/output/merged.tree.final.nwk"
tax_file = directory+"/output/merged.taxonomy.final.tsv"
count_table = directory+"/output/merged.counts.final.biom"
# Export Classifications
os.system(' '.join([
    "qiime tools export",
    directory+"/output/merged.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/merged.table.final.qza",
    "--output-dir "+directory+"/output/"
]))
# Export SV Sequences
os.system(' '.join([
    "qiime tools export",
    directory+"/output/merged.rep.seqs.final.qza",
    "--output-dir "+directory+"/output/"
]))
if domain == "bacteria":
    # Export Tree
    os.system(' '.join([
        "qiime tools export",
        directory+"/output/merged.rep.seqs.tree.final.qza",
        "--output-dir "+directory+"/output/"
    ]))
    %mv $directory/output/tree.nwk $tree_file
# 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
Explanation: Step 6: Prepare Data for Import to Phyloseq
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)
#### 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/merged.taxonomy.final.tsv", sep="\t")
%R -i tax_file
%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")
%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/merged.counts.final.biom",
    "-o",
    directory+"/output/merged.counts.final.tsv",
    "--to-tsv"
]))
# The biom converter adds a stupid line that messes with the table formatting
fix_biom_conversion(directory+"/output/merged.counts.final.tsv")
# Finally import
count_table = pd.read_csv(directory+"/output/merged.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)
if domain == "bacteria":
    # Import Phylogenetic Tree
    tree_file = directory+"/output/merged.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_merged.final.rds"
%R -i output
%R saveRDS(p_final, file = output)
# Confirm Output
%R print(p_final)
Explanation: Step 7: Import into Phyloseq
End of explanation
# Remove Files
if domain == "bacteria":
    %rm $directory/output/*tree.unrooted.qza 
    %rm $directory/output/*aligned.masked.qza 
    
%rm $directory/output/*.biom 
%rm $directory/output/taxonomy.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/
# 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 directories 'final' and 'intermediate_files'")
Explanation: Step 8: Clean-up Intermediate Files and Final Outputs
End of explanation
<END_TASK>