<commit_before># <NAME>

# This pipeline processes metagenomic reads, assembles, and annotates them to enable metatranscriptomic mapping

# Table of Contents
# 1. Combine and trim metagenomic reads
# 2. 
# 3. 
# 4. 
# 5. 
# 6. 
# 7. 
# 8. 
# 9. 
# 10. 

#-------------------------------------------------------------------------------------------------------------------------#

# 1. Process raw metagenomic read data

# Combine files from the same lane and read
cat *L001_R1* > metagenome.read1.fastq
cat *L001_R2* > metagenome.read2.fastq

# Interleave the paired reads
source /home/mljenior/bin/khmerEnv/bin/activate
python /home/mljenior/bin/khmerEnv/bin/interleave_reads.py metagenome.read1.fastq metagenome.read1.fastq metagenome.paired.fastq

# Quality trim with sickle (ends and ambiguous bases)
#/mnt/EXT/Schloss-data/kiverson/sickle-master/sickle pe -c metagenome.paired.fastq -M metagenome.paired.trimmed.fastq -t sanger

# Return some information about the quality trimmed, paired-end fastq
python /mnt/EXT/Schloss-data/seq_stats.py metagenome.paired.trimmed.fastq > metagenome.paired.trimmed.fastq.stats.txt


# In parallel, convert directly to fasta format to map to contigs later
awk '{print ">" substr($0,2);getline;print;getline;getline}' metagenome.paired.trimmed.fastq > metagenome.paired.trimmed.fasta

#-------------------------------------------------------------------------------------------------------------------------#

# 2. Normalization with khmer

# 2-pass Digital normalization
source /home/mljenior/bin/khmerEnv/bin/activate
python /home/mljenior/bin/khmerEnv/bin/normalize-by-median.py -p -k 20 -C 20 -N 4 -x 8e9 --savetable normC20k20.kh metagenome.paired.trimmed.fastq
python /home/mljenior/bin/khmerEnv/bin/filter-abund.py -V normC20k20.kh metagenome.paired.trimmed.fastq.keep
python /home/mljenior/bin/khmerEnv/bin/extract-paired-reads.py metagenome.paired.trimmed.fastq.keep.abundfilt
python /home/mljenior/bin/khmerEnv/bin/normalize-by-median.py -C 5 -k 20 -N 4 -x 8e9 --savetable normC5k20.kh -p metagenome.paired.trimmed.fastq.keep.abundfilt.pe

# Convert the normalized file to a renamed fasta
awk '{print ">" substr($0,2);getline;print;getline;getline}' metagenome.paired.trimmed.fastq.keep.abundfilt.pe.keep > metagenome.paired.trimmed.normalized.fasta

# Return some more information about the results of digital normalization
python /mnt/EXT/Schloss-data/seq_stats.py metagenome.paired.trimmed.normalized.fasta > metagenome.paired.trimmed.normalized.fasta.stats.txt

#-------------------------------------------------------------------------------------------------------------------------#

# 3. Clean up

# Delete excess files created during normalization steps
rm *.keep *.abundfilt *.pe *.se *.kh metagenome.paired.trimmed.fasta

# Move the 2 files you will need later
mv *.fasta ../assembly

# Zip files that are not going to be used downstream
qzip *fastq

#-------------------------------------------------------------------------------------------------------------------------#

# 4. Assembly - In new directory

# Assemble using Megahit
python /home/mljenior/bin/megahit/megahit -m 45e9 -l 251 -r metagenome.paired.trimmed.normalized.fasta --cpu-only -o metagenome.paired.trimmed.normalized.megahit

# Remove contigs shorter than the read size
python /share/scratch/bin/removeShortContigs.py metagenome.paired.trimmed.normalized.megahit/final.contigs.fa metagenome.paired.trimmed.normalized.megahit/final.contigs.250.fa 250

# Get some information on the assembled contigs
python /mnt/EXT/Schloss-data/seq_stats.py metagenome.paired.trimmed.normalized.megahit/final.contigs.250.fa > metagenome.paired.trimmed.normalized.megahit/final.contigs.250.fa.stats.txt

#-------------------------------------------------------------------------------------------------------------------------#

# 5. QC of assembly

# Check how well the reads (uniques from first section) map back to the newly assembled contigs
/home/mljenior/bin/bowtie/bowtie-build metagenome.paired.trimmed.normalized.megahit/final.contigs.250.fa megahit_database
/home/mljenior/bin/bowtie/bowtie megahit_database -f metagenome.paired.trimmed.unique.fasta -p 4 -S aligned_reads.megahit.sam

# The # of reads mapping to contigs will print to the screen

#-------------------------------------------------------------------------------------------------------------------------#

# 6. Call genes and return fasta files of sequences

# Find genes using mga
/share/scratch/bin/mga final.contigs.250.fa -m > final.contigs.250.mga.out

# script to get genes from mga.out and write them to amino acid and nucleotide fastas
python /share/scratch/bin/getgenesbypos.py final.contigs.250.mga.out final.contigs.250.fa genes.faa genes.fna

# remove short genes
python /share/scratch/bin/removeShortContigs.py genes.fna genes.250.fna 250
python /share/scratch/bin/removeShortContigs.py genes.faa genes.80.faa 80

# Get some information on the called gene sequences
python /mnt/EXT/Schloss-data/seq_stats.py genes.250.fna > genes.250.fna.stats.txt

#-------------------------------------------------------------------------------------------------------------------------#

# 7. Make KEGG BLAST database








#-------------------------------------------------------------------------------------------------------------------------#


# 8. Annotate called genes

# Blast genes
blastn genes.250.fna






# Reassociate BLAST results with gene sequences




# Create Bowtie database of annotated genes
/home/mljenior/bin/bowtie/bowtie-build annotated.contigs.250.fa metagenome_megahit_database


#-------------------------------------------------------------------------------------------------------------------------#

9. Process raw metatransciptomic reads

# Combine files from the same lane and read
cat *L001_R1* > metatransciptome.read1.fastq
cat *L001_R2* > metatransciptome.read2.fastq

# Interleave the paired reads
source /home/mljenior/bin/khmerEnv/bin/activate
python /home/mljenior/bin/khmerEnv/bin/interleave_reads.py metatransciptome.read1.fastq metatransciptome.read1.fastq metatransciptome.paired.fastq

# Quality trim with sickle (ends and ambiguous bases)
#/mnt/EXT/Schloss-data/kiverson/sickle-master/sickle pe -c metatransciptome.paired.fastq -M metatransciptome.paired.trimmed.fastq -t sanger

# Convert the fastq to a fasta
awk '{print ">" substr($0,2);getline;print;getline;getline}' metatransciptome.paired.trimmed.fastq > metatransciptome.paired.trimmed.fasta

# Get some information on the transcript file
python /mnt/EXT/Schloss-data/seq_stats.py metatransciptome.paired.trimmed.fasta > metatransciptome.paired.trimmed.fasta.stats.txt

# Move the processed reads to a new directory
mv metatransciptome.paired.trimmed.fasta ../transcript_mapping

# Clean up
gzip *fastq

#-------------------------------------------------------------------------------------------------------------------------#

# 10. Map transcripts to genes
/home/mljenior/bin/bowtie/bowtie metagenome_megahit_database -f metatransciptome.paired.trimmed.fasta -p 4 -S aligned_reads.megahit.sam



# Convert the sam to bam and get rid of the old sam file.
samtools view -bS alignedreads.sam > alignedreads.bam
rm alignedreads.sam


# Sort the bamfile
samtools sort alignedreads.bam alignedreads.sorted


#create an index. This will make searching for things in the file faster
samtools index alignedreads.sorted.bam


# use the index to give you stats about your bam in human readable format. This is where you can see how many reads mapped to each contig
samtools idxstats alignedreads.sorted.bam > $i.mapped2contigs.txt



<commit_msg>changed assembly to take paired and single end reads<commit_after># <NAME>

# This pipeline processes metagenomic reads, assembles, and annotates them to enable metatranscriptomic mapping

# Table of Contents
# 1. Combine and trim metagenomic reads
# 2. 
# 3. 
# 4. 
# 5. 
# 6. 
# 7. 
# 8. 
# 9. 
# 10. 

#-------------------------------------------------------------------------------------------------------------------------#

# 1. Process raw metagenomic read data

# Combine files from the same lane and read
cat *L001_R1* > metagenome.read1.fastq
cat *L001_R2* > metagenome.read2.fastq

# Interleave the paired reads
source /home/mljenior/bin/khmerEnv/bin/activate
python /home/mljenior/bin/khmerEnv/bin/interleave_reads.py metagenome.read1.fastq metagenome.read1.fastq metagenome.paired.fastq

# Quality trim with sickle (ends and ambiguous bases)
/mnt/EXT/Schloss-data/kiverson/sickle-master/sickle pe -c metagenome.paired.fastq -M metagenome.paired.trimmed.fastq -t sanger

# Return some information about the quality trimmed, paired-end fastq
python /mnt/EXT/Schloss-data/seq_stats.py metagenome.paired.trimmed.fastq > metagenome.paired.trimmed.fastq.stats.txt


# In parallel, convert directly to fasta format to map to contigs later and move it to the correct directory
awk '{print ">" substr($0,2);getline;print;getline;getline}' metagenome.paired.trimmed.fastq > metagenome.paired.trimmed.fasta
mv metagenome.paired.trimmed.fasta ../assembly/bowtie

#-------------------------------------------------------------------------------------------------------------------------#

# 2. Normalization with khmer

# 2-pass Digital normalization
source /home/mljenior/bin/khmerEnv/bin/activate
python /home/mljenior/bin/khmerEnv/bin/normalize-by-median.py -p -k 20 -C 20 -N 4 -x 8e9 --savetable normC20k20.kh metagenome.paired.trimmed.fastq
python /home/mljenior/bin/khmerEnv/bin/filter-abund.py -V normC20k20.kh metagenome.paired.trimmed.fastq.keep
python /home/mljenior/bin/khmerEnv/bin/extract-paired-reads.py metagenome.paired.trimmed.fastq.keep.abundfilt
python /home/mljenior/bin/khmerEnv/bin/normalize-by-median.py -C 5 -k 20 -N 4 -x 8e9 --savetable normC5k20.kh -p metagenome.paired.trimmed.fastq.keep.abundfilt.pe

# Convert the normalized file to a renamed fasta
awk '{print ">" substr($0,2);getline;print;getline;getline}' metagenome.paired.trimmed.fastq.keep.abundfilt.pe.keep > metagenome.trimmed.normalized.paired.fasta

# Return some more information about the results of digital normalization
python /mnt/EXT/Schloss-data/seq_stats.py metagenome.trimmed.normalized.paired.fasta > metagenome.trimmed.normalized.paired.fasta.stats.txt

#-------------------------------------------------------------------------------------------------------------------------#

# 3. Clean up

# Move the 2 files you will need later
mv *.fasta ../assembly

# Delete excess files created during normalization steps
rm *.keep *.abundfilt *.pe *.se *.kh *.paired.trimmed.fasta

# Zip files that are not going to be used downstream
gzip *fastq

#-------------------------------------------------------------------------------------------------------------------------#

# 4. Assembly - In new directory

# Assemble using Megahit
python /home/mljenior/bin/megahit/megahit -m 45e9 -l 251 --12 metagenome.trimmed.normalized.paired.fasta -r metagenome.trimmed.normalized.single.fasta --cpu-only --k-max 127 -o metagenome.paired.trimmed.normalized.megahit

# Remove contigs shorter than the read size
python /share/scratch/bin/removeShortContigs.py metagenome.paired.trimmed.normalized.megahit/final.contigs.fa metagenome.paired.trimmed.normalized.megahit/final.contigs.250.fa 250

# Get some information on the assembled contigs
python /mnt/EXT/Schloss-data/seq_stats.py metagenome.paired.trimmed.normalized.megahit/final.contigs.250.fa > metagenome.paired.trimmed.normalized.megahit/final.contigs.250.fa.stats.txt

#-------------------------------------------------------------------------------------------------------------------------#

# 5. QC of assembly

# Check how well the reads (uniques from first section) map back to the newly assembled contigs
/home/mljenior/bin/bowtie/bowtie-build metagenome.paired.trimmed.normalized.megahit/final.contigs.250.fa metagenome_megahit_database
/home/mljenior/bin/bowtie/bowtie metagenome_megahit_database -f metagenome.paired.trimmed.fasta -p 4 -S aligned_reads.megahit.sam

# The # and % of reads mapping to contigs will print to the screen

#-------------------------------------------------------------------------------------------------------------------------#

# 6. Call genes and return fasta files of sequences

# Find genes using mga
/share/scratch/bin/mga final.contigs.250.fa -m > final.contigs.250.mga.out

# script to get genes from mga.out and write them to amino acid and nucleotide fastas
python /share/scratch/bin/getgenesbypos.py final.contigs.250.mga.out final.contigs.250.fa genes.faa genes.fna

# remove short genes
python /share/scratch/bin/removeShortContigs.py genes.fna genes.250.fna 250
python /share/scratch/bin/removeShortContigs.py genes.faa genes.80.faa 80

# Remove duplicate sequences
/mnt/EXT/Schloss-data/bin/fastx_collapser -i genes.250.fna -o genes.250.unique.fna
/mnt/EXT/Schloss-data/bin/fastx_collapser -i genes.80.faa -o genes.80.unique.faa

# Get some information on the called gene sequences
python /mnt/EXT/Schloss-data/seq_stats.py genes.250.unique.fna > genes.250.unique.fna.stats.txt
python /mnt/EXT/Schloss-data/seq_stats.py genes.80.unique.faa > genes.80.unique.faa.stats.txt

#-------------------------------------------------------------------------------------------------------------------------#

# 7. Make KEGG BLAST database

makeblastdb -in /kegg/genes/fasta/genes.pep -dbtype 'prot' -out kegg_prot_blast.db

#-------------------------------------------------------------------------------------------------------------------------#

# 8. Annotate called genes

# Blast genes (all 6 frames)
blastp -query genes.80.unique.faa -db kegg_prot_blast.db -num_threads 16 -evalue 1e-10 -best_hit_score_edge 0.05 -best_hit_overhang 0.25 -outfmt 6 -max_target_seqs 1 -out metagenome.protVprot.out

# Associate BLAST results with gene sequences
python /home/mljenior/scripts/annotate_fasta.py genes.80.unique.faa metagenome.protVprot.out > annotated.genes.80.unique.faa

# Create BLAST database of annotated genes
makeblastdb -in ../assembly/annotated.genes.80.unique.faa -dbtype 'prot' -out annotated_metagenome_blast.db

#-------------------------------------------------------------------------------------------------------------------------#

9. Process raw metatransciptomic reads

# Combine files from the same lane and read
cat *L001_R1* > metatransciptome.read1.fastq
cat *L001_R2* > metatransciptome.read2.fastq

# Interleave the paired reads
source /home/mljenior/bin/khmerEnv/bin/activate
python /home/mljenior/bin/khmerEnv/bin/interleave_reads.py metatransciptome.read1.fastq metatransciptome.read1.fastq metatransciptome.paired.fastq

# Quality trim with sickle (ends and ambiguous bases)
#/mnt/EXT/Schloss-data/kiverson/sickle-master/sickle pe -c metatransciptome.paired.fastq -M metatransciptome.paired.trimmed.fastq -t sanger

# Convert the fastq to a fasta
awk '{print ">" substr($0,2);getline;print;getline;getline}' metatransciptome.paired.trimmed.fastq > metatransciptome.paired.trimmed.fasta

# Get some information on the transcript file
python /mnt/EXT/Schloss-data/seq_stats.py metatransciptome.paired.trimmed.fasta > metatransciptome.paired.trimmed.fasta.stats.txt

# Move the processed reads to a new directory
mv metatransciptome.paired.trimmed.fasta ../transcript_mapping

# Clean up
gzip *fastq

#-------------------------------------------------------------------------------------------------------------------------#

# 10. Map transcripts to genes

tblastn -query metatransciptome.paired.trimmed.fasta -db annotated_metagenome_blast.db -num_threads 16 -evalue 1e-10 -best_hit_score_edge 0.05 -best_hit_overhang 0.25 -outfmt 6 -max_target_seqs 1 -out mapped_transcriptsVmetagenome.out


Write a script to create a table from the BLAST results

{'fraction_non_alphanumeric': 0.23743511836941386, 'fraction_numerical': 0.020129130269654388, 'mean_word_length': 2.230879345603272, 'pattern_counts': {'":': 0, '<': 5, '<?xml version=': 0, '>': 41, '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}}
{'dir': 'git-commits-cleaned', 'id': '5481071', 'n_tokens_mistral': 5225, 'n_tokens_neox': 4742, 'n_words': 1158}