<jupyter_start><jupyter_text>Notebook for processing data as described in 'RNA viral communities are structured by host plant phylogeny in oak and conifer leaves'<NAME> Horst  QC reads with trimmomaticTrim the reads with trimmomatic, following Roux et al., 2017minimum base quality threshold of 30 evaluated on sliding windows of 4 bases, and minimum read length of 50As a rule of thumb newer libraries will use TruSeq3, but this really depends on your service provider.http://www.usadellab.org/cms/?page=trimmomatic<jupyter_code># loading modules
module load java
module load trimmomatic

# use trimmomatic to trim adapters, make sure you have the correct adapterfile
trimmomatic PE -threads 8 -phred33 $1 ${1%%_1*}_2.fastq.gz \
../trimmed/${1%%_1*}_1_trimmed.fq.gz unpaired/${1%%_1*}_1_Unpaired.fq.gz \
../trimmed/${1%%_1*}_2_trimmed.fq.gz unpaired/${1%%_1*}_2_Unpaired.fq.gz \
ILLUMINACLIP:/adapters/TruSeq3-PE.fa:2:30:10 \
SLIDINGWINDOW:4:30 MINLEN:50<jupyter_output><empty_output><jupyter_text>Remove PhiX174 with bbdukPhix174 is a phage sequence that is used to aid in sequencing. This needs to be removed from our sequencing datahttps://jgi.doe.gov/data-and-tools/bbtools/bb-tools-user-guide/bbduk-guide/<jupyter_code># load module
module load java
module load bbmap

# use bbduk to remove phix contamination
bbduk.sh in1=$1 in2=${1%%_1*}_2_trimmed.fq.gz \
out1=../remove_phix/${1%%_1*}_1_trimmed_bbduk.fq.gz out2=../remove_phix/${1%%_1*}_2_trimmed_bbduk.fq.gz \
ref=/bbduk/phix_genome.fa k=31 \
hdist=1 stats=../remove_phix/${1%%}_stats.txt<jupyter_output><empty_output><jupyter_text>Assembly with MEGAHITAssemble clean reads using Megahit. Use standard settings, which menas the minimum length of a contig will be 200 bphttps://github.com/voutcn/megahit<jupyter_code># loading modules
module load megahit

# Run megahit using standard settings
megahit -1 $1 -2 ${1%%1_R1*}1_R2_trimmed_bbduk.fq.gz \
-r ../unpaired/${1%%1_R1*}_R1_unpaired.fq.gz,../unpaired/${1%%1_R1*}_R2_unpaired.fq.gz \
-o ../assemblies_persample/${1} \
--out-prefix ${1}<jupyter_output><empty_output><jupyter_text>Search for RdRp sequences in assemblies data- Following Starr et al., 2019- Use PFAM numbers Mononeg_RNA_pol [PF00946], RdRP_5 [PF07925], Flavi_NS5 [PF00972], Bunya_RdRp [PF04196], Mitovir_RNA_pol [PF05919], RdRP_1 [PF00680], RdRP_2 [PF00978], RdRP_3 [PF00998], RdRP_4 [PF02123], Viral_RdRp_C [PF17501], and Birna_RdRp [PF04197]. - Download Seeds of the HMM alignments and use hmmer to make HMM profiles for each- Predict proteins with Prodigal on assemblies- Use Hmmer to search predicted genes for RdRps<jupyter_code># First rename the assemblies so results are easier to look at

# Using bbmap
module load bbmap

# rename all contigs in each fasta file after the fasta file
for f in *.fa
do 
echo $f
rename.sh in=$f out=${f%%.fa*}_renamed.fa prefix=${f%%.contigs.fa*}
done<jupyter_output><empty_output><jupyter_text>PFAM numbers for RdRps- Read more about RdRps here: http://pfam.xfam.org/family/PF00680- Download the alignment in selex format from the PFAM website- Do not use PFAMs for RVT_1 and RVT_2, because retrotransposons have those too.<jupyter_code># build Hmm profiles from the selex files
# load module
module load hmmer

# loop through selex files
for f in *.txt
do
# Use hmmbuild to make hmm profile for each RdRp
hmmbuild ${f%%_seed*}.hmm $f
done 

# Concatenate all hmmprofiles so we only have to do one search
cat *.hmm >> all_RdRp.hmm<jupyter_output><empty_output><jupyter_text>Prodigal for gene predictions- USe the anonymous mode (or meta mode, is the same)- Read about why meta mode: https://www.biocode.ltd/catalog-tooldata/prodigal- Do HMM search on these predictions- Cluster the found predicted RdRp protein sequences using USEARCH<jupyter_code># Load modules
module load prodigal
module load hmmer

# Loop through fasta files (contigs)
for f in *.fa
do

# predict genes with prodigal for each file
prodigal -i $f -a ../prodigal/${f%%.contigs.fa*}_proteins.faa -p meta -q
done

# Do an HMM search on these prodigal gene predictions
# Loop through the gene prediction files (proteins)
for f in *.faa
do
echo $f
# Use Hmmsearch to look for RdRp sequences
hmmsearch -E 0.00001 -o ../hmm_search/out/${f%%_proteins.faa*}.out \
--tblout ../hmm_search/${f%%_proteins.faa*}_tabular_eval.out  \
../hmms/all_RdRp.hmm $f 
done

# Cluster RdRp sequences with USEARCH
../../programs/usearch -cluster_fast \
./prodigal/predicted_rdrp.faa -id 0.99 -centroids \
./prodigal/predicted_rdrp_clust.faa<jupyter_output><empty_output><jupyter_text>Run VIBRANT on all contigs - Changed VIBRANT source code so it will take contigs of 100 bp with one orf- changed /share/emersonlab/annie/miniconda3/envs/vibrant/bin/VIBRANT_run.py- Try on all contigs<jupyter_code>## Changed line 139 in VIBRANT_run.py as follows:
if int(orf_low) < 1 or int(lim_low) < 100:
        logging.info("VIBRANT error: minimum ORFs is 4 and minimum sequence length is 1kb. These variables can only increase. Exiting." + "\n")
        exit()
    command = ' -f ' + str(format) + ' -d ' + str(databases) + ' -m ' + str(files) + ' -l ' + str(lim_low) + ' -o ' + str(orf_low) + str(virome)
    subprocess.run("echo '' > " + str(out_folder)+base + "_four-orf-count.txt", shell=True)

# Ran this adapted version of vibrant as follows:
VIBRANT_run.py -i all_contigs.fa -folder vibrant -t 5 -l 100 -o 1<jupyter_output><empty_output><jupyter_text>Take all contigs with predicted viral protein- VIBRANT output that it predicted to be viral- RdRP matches to proteins- Concatenate these files- remove duplicates- predict proteins- Run blastP to complete nr database- Run blastN to complete nt database Use code in the 210120_concat_alloutputs notebook to concatenate:- Hmmr output- Vibrant output- BlastP output- BlastN output<jupyter_code># remove duplicate contigs using bbmap
dedupe.sh in=210112_all_virallike_contigs.fa \
out=210112_all_virallike_contigs_dedup.fa ac=f requirematchingnames

# predict proteins
module load prodigal
prodigal -i 210112_all_virallike_contigs_dedup.fa \
-a 210112_all_virallike_contigs_dedup.faa -p meta -q

# remove spaces from protein fasta headers
sed '/^>/ s/ .*//' 210112_all_virallike_contigs_dedup.faa >> 210112_all_virallike_contigs_dedup_nospaces.faa

# Run BlastP 
blastp -db /share/genomes/ncbi/nr -query 210112_all_virallike_contigs_dedup_nospaces.faa \
-out ../blastp_out/210112_all_virallike_contigs_blastP.csv \
-outfmt 10 -max_target_seqs 5 -num_threads 15 

# Run BlastN
blastn -db /share/genomes/ncbi/nt -query 210112_all_virallike_contigs_dedup.fa \
-out ../210112_all_virallike_contigs_blastN.csv \
-outfmt 10 -max_target_seqs 5 -num_threads 15<jupyter_output><empty_output><jupyter_text>Dereplicate the contigs using dRephttps://drep.readthedocs.io/en/latest/<jupyter_code># dereplicate with dRep
dRep dereplicate ./drep --S_algorithm ANImf -sa 0.95 -nc 0.85 -g ./contigs/*.fa -p 6<jupyter_output><empty_output><jupyter_text>Map back the clean reads to contigs- Contigs were manually curated- Making sure that there were no retroviral elements in there, cause these are likely to be retrotransposons- Use Bowtie2 to map back- Use bbmap and samtools to index the bam files<jupyter_code>#  make reference file with bowtie
module load bowtie2
bowtie2-build --threads 4 viral_like_contigs.fa viral_like_contigs

# map back clean reads with bowtie
module load bowtie2

# For each read file do bowtie2, using sensitive preset
for f in *_R1_*.fq.gz
do
bowtie2 --threads 3 \
-x viral_like_contigs \
-1 $f \
-2 ${f%%_R1_*}_R2_trimmed_bbduk.fq.gz \
-S  ${f%%_R1_*}.sam --no-unal --sensitive

# index and convert sam to bam
module load bbmap
module load samtools

for f in *.sam
do echo $f
samtools view -F 4 -bS $f | samtools sort > ${f%%.sam*}_sortedIndexed.bam
done

for f in *.bam ; do samtools index $f ; done;<jupyter_output><empty_output><jupyter_text>Use coverM to make a coverage table- Creates abundance tables- https://github.com/wwood/CoverM- Using 75% coverage<jupyter_code># Run CoverM
coverm contig -m mean --min-covered-fraction 0.75 -b *.bam > coverage_table_arboretum_75_mean.csv
coverm contig -m count --min-covered-fraction 0.75 -b *.bam > coverage_table_arboretum_75_count.csv<jupyter_output><empty_output><jupyter_text>Make phylogenetic trees of the RdRp sequences- Dereplicate RdRp protein sequences using CDhit- Phylogenetic alignments using Mafft- remove spurious seqs using trimal- infer phylogenetic trees using FastTree<jupyter_code># remove the spaces from own RdRp sequences
sed '/^>/ s/ .*//' clustered_predict_rdrp.faa >> \
210329_clustered_predict_rdrp_nospace.faa

# remove predicted RdRp sequences shorter than 60 AA, cause they will interfere with alignment
# Manually, removed 7 seqs

# cluster the RdRps from refseq with USEARCH
# 99% AAI (Starr et al., 2019), 631 clusters
../../programs/usearch -cluster_fast 210329_RdRp_refseq_mancur.faa \
-id 0.99 -centroids 210329_RdRp_refseq_clust.fasta 


# align with MAFFT (Shi2016)
# first make file with all RdRps
cat 210329_clustered_predict_rdrp_nospace.faa 210329_RdRp_refseq_clust.fasta >> 210329_rdrp_own_refseq.faa

# Use auto setting, leave gappyregions
module load mafft
mafft --thread 5 \
--adjustdirection \
--auto \
--reorder \
210329_rdrp_own_refseq.faa > 210329_rdrp_own_refseq_mafftauto.faa

# Use trimal to remove spurious sequences
trimal -in 210329_rdrp_own_refseq_mafftauto.faa -out 210329_rdrp_own_refseq_mafftauto_trimal.faa -gappyout

# Use FastTree to make phylogenetic trees
FastTreeMP < 210329_rdrp_own_refseq_mafftauto_trimal.faa > fasttree_tree.nwk \
-wag -log 210516_fasttree.log<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.09772325605572305, 'fraction_numerical': 0.040128911529264996, 'mean_word_length': 4.873015873015873, 'pattern_counts': {'":': 0, '<': 53, '<?xml version=': 0, '>': 66, 'https://': 5, 'lorem ipsum': 0, 'www.': 2, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '12', 'dir': 'jupyter-structured-clean-dedup', 'id': '603387', 'license': 'MIT', 'n_tokens_mistral': 3875, 'n_tokens_neox': 3467, 'n_words': 1018, 'path': 'Data_processing_arboretum.ipynb', 'repo_name': 'AnneliektH/arboretum', 'size': '9625'}