Monkeypox has been termed a global health emergency. For researchers interested in tracking the evolution of the monkeypox virus genome, it can be difficult to get sufficient coverage. The xGen Monkeypox Virus Amplicon Panel provides a single-tube, two-step PCR amplification workflow with primer sets designed to create amplicons across the entire genome.*
If coverage is too low, you may need to sequence more reads to ensure uniform coverage. For many assembly algorithms, 30x coverage or higher is recommended.
Query Seq-id Start of alignment in query End of alignment in query Start of alignment in subject End of alignment in subject Expect value Alignment length Percentage of identical matches Subject Seq-id Subject Title Query Coverage Per Subject Query Coverage Per HSP Subject accession Subject sequence length Query sequence length
-
Using bacto pipeline to get trimmed reads
cp /home/jhuang/Tools/bacto/bacto-0.1.json . cp /home/jhuang/Tools/bacto/cluster.json . cp /home/jhuang/Tools/bacto/Snakefile . ln -s /home/jhuang/Tools/bacto/local . ln -s /home/jhuang/Tools/bacto/db . ln -s /home/jhuang/Tools/bacto/envs . mkdir raw_data; cd raw_data ln -s ../raw_data_downloads/Affe30_S1_L001_R1_001.fastq.gz Affe30_R1.fastq.gz ln -s ../raw_data_downloads/Affe30_S1_L001_R2_001.fastq.gz Affe30_R2.fastq.gz ln -s ../raw_data_downloads/Affe31_S2_L001_R1_001.fastq.gz Affe31_R1.fastq.gz ln -s ../raw_data_downloads/Affe31_S2_L001_R2_001.fastq.gz Affe31_R2.fastq.gz (bengal3_ac3) /home/jhuang/miniconda3/envs/snakemake_4_3_1/bin/snakemake --printshellcmds -
Filtering low complexity
fastp -i Affe30_trimmed_P_1.fastq -I Affe30_trimmed_P_2.fastq -o Affe30_trimmed_R1.fastq -O Affe30_trimmed_R2.fastq --low_complexity_filter --complexity_threshold 30 fastp -i Affe31_trimmed_P_1.fastq -I Affe31_trimmed_P_2.fastq -o Affe31_trimmed_R1.fastq -O Affe31_trimmed_R2.fastq --low_complexity_filter --complexity_threshold 30 -
Using vrap to assembly and annotate the contigs, the spades-step was replaced with idba of DAMIAN
#--host /home/jhuang/REFs/genome.fa -n /mnt/nvme0n1p1/blast/nt -a /mnt/nvme0n1p1/blast/nr vrap/vrap.py -1 Affe30_trimmed_R1.fastq.gz -2 Affe30_trimmed_R2.fastq.gz -o Affe30_trimmed_vrap_out -t 40 -l 100 vrap/vrap.py -1 Affe31_trimmed_R1.fastq.gz -2 Affe31_trimmed_R2.fastq.gz -o Affe31_trimmed_vrap_out -t 40 -l 100 #--> ERROR in spades-assembly, we usding idba from DAMIAN assembly, copy the assembly to spades. IT WORKS! damian.rb --host human3 --type dna -1 /home/jhuang/DATA/Data_Susanne_MPox/Affe31_trimmed_R1.fastq.gz -2 /home/jhuang/DATA/Data_Susanne_MPox/Affe31_trimmed_R2.fastq.gz --sample Affe31_megablast --blastn never --blastp never --min_contiglength 100 --threads 56 --force damian.rb --host human3 --type dna -1 /home/jhuang/DATA/Data_Susanne_MPox/Affe30_trimmed_R1.fastq.gz -2 /home/jhuang/DATA/Data_Susanne_MPox/Affe30_trimmed_R2.fastq.gz --sample Affe30_megablast --blastn never --blastp never --min_contiglength 100 --threads 56 --force damian.rb --host human3 --type dna -1 /home/jhuang/DATA/Data_Susanne_MPox/Affe31_trimmed_R1.fastq.gz -2 /home/jhuang/DATA/Data_Susanne_MPox/Affe31_trimmed_R2.fastq.gz --sample Affe31_blastn --blastn progressive --blastp never --min_contiglength 100 --threads 56 --force damian.rb --host human3 --type dna -1 /home/jhuang/DATA/Data_Susanne_MPox/Affe30_trimmed_R1.fastq.gz -2 /home/jhuang/DATA/Data_Susanne_MPox/Affe30_trimmed_R2.fastq.gz --sample Affe30_blastn --blastn progressive --blastp never --min_contiglength 100 --threads 56 --force cp ~/rtpd_files/Affe30_megablast/idba_ud_assembly/contig.fa contigs.fasta cp ~/rtpd_files/Affe31_megablast/idba_ud_assembly/contig.fa contigs.fasta vrap/vrap.py -1 Affe30_trimmed_R1.fastq.gz -2 Affe30_trimmed_R2.fastq.gz -o Affe30_trimmed_vrap_out -t 40 -l 100 vrap/vrap.py -1 Affe31_trimmed_R1.fastq.gz -2 Affe31_trimmed_R2.fastq.gz -o Affe31_trimmed_vrap_out -t 40 -l 100 -
Step-by-Step Process for read and contig mapping profiling plots
a. mapping the contig on the reference JX878414 bowtie2-build JX878414.1.fasta JX878414.1_index bowtie2 -f -x JX878414.1_index -U Affe30_trimmed_vrap_out/spades/contigs.fasta -S Affe30_contigs_aligned.sam samtools view -S -b Affe30_contigs_aligned.sam > Affe30_contigs_aligned.bam samtools sort Affe30_contigs_aligned.bam -o Affe30_contigs_aligned_sorted.bam samtools index Affe30_contigs_aligned_sorted.bam bowtie2 -f -x JX878414.1_index -U Affe31_trimmed_vrap_out/spades/contigs.fasta -S Affe31_contigs_aligned.sam samtools view -S -b Affe31_contigs_aligned.sam > Affe31_contigs_aligned.bam samtools sort Affe31_contigs_aligned.bam -o Affe31_contigs_aligned_sorted.bam samtools index Affe31_contigs_aligned_sorted.bam 62 reads; of these: 62 (100.00%) were unpaired; of these: 1 (1.61%) aligned 0 times 61 (98.39%) aligned exactly 1 time 0 (0.00%) aligned >1 times 98.39% overall alignment rate 17 reads; of these: 17 (100.00%) were unpaired; of these: 1 (5.88%) aligned 0 times 16 (94.12%) aligned exactly 1 time 0 (0.00%) aligned >1 times 94.12% overall alignment rate b. candidate methods for mapping the contigs on a reference #1. Minimap2 – Efficient for aligning large contigs or entire genomes. minimap2 -ax asm5 JX878414.1.fasta contigs_Affe30.fasta > output.sam #2. BLAST – More sensitive but slower; good for divergent sequences. blastn -query contigs_Affe30.fasta -db JX878414.1.fasta -out output.blast -outfmt 6 #3. LAST – Suitable for more distant relationships and large genomic changes. lastdb my_reference_db JX878414.1.fasta lastal my_reference_db contigs_Affe30.fasta > output.maf #4. MUMmer (nucmer) – Best for comparing large contigs or entire genomes, especially with structural variations. nucmer JX878414.1.fasta contigs_Affe30.fasta -p output c. Index the Reference Genome (JX878414.1.fasta): You'll need to index your reference genome to align the contigs and reads against it. samtools faidx JX878414.1.fasta d. Generate Coverage Profile for Reads (from Fastq): Align the trimmed fastq reads (Affe31_trimmed_R1.fastq.gz and Affe31_trimmed_R2.fastq.gz) to the reference genome using a mapper like BWA or Bowtie2. bwa index JX878414.1.fasta bwa mem JX878414.1.fasta Affe31_trimmed_R1.fastq.gz Affe31_trimmed_R2.fastq.gz > Affe31_reads_aligned.sam Convert the SAM file to BAM and sort it: samtools view -Sb Affe31_reads_aligned.sam | samtools sort -o Affe31_reads_aligned_sorted.bam samtools index Affe31_reads_aligned_sorted.bam #6743 + 0 in total (QC-passed reads + QC-failed reads) #0 + 0 secondary #9 + 0 supplementary #0 + 0 duplicates #5312 + 0 mapped (78.78% : N/A) #6734 + 0 paired in sequencing #3367 + 0 read1 #3367 + 0 read2 #4822 + 0 properly paired (71.61% : N/A) #4844 + 0 with itself and mate mapped #459 + 0 singletons (6.82% : N/A) #0 + 0 with mate mapped to a different chr #0 + 0 with mate mapped to a different chr (mapQ>=5) e. Generate Coverage Profile for Contigs: The file Affe31_output_sorted.bam is already aligned against the reference genome. You can directly use it to visualize contig coverage. Ensure the BAM file is indexed: samtools index Affe31_contigs_aligned_sorted.bam f. Generate Coverage Tracks: Use BamCoverage to generate coverage files (in bigWig format) for both the reads and contigs. For reads coverage: bamCoverage -b Affe31_reads_aligned_sorted.bam -o Affe31_reads_coverage.bw For contigs coverage: bamCoverage -b Affe31_contigs_aligned_sorted.bam -o Affe31_contigs_coverage.bw g. Generate Coverage Profile for Reads (from Fastq): Align the trimmed fastq reads (Affe30_trimmed_R1.fastq.gz and Affe30_trimmed_R2.fastq.gz) to the reference genome using a mapper like BWA or Bowtie2. bwa index JX878414.1.fasta bwa mem JX878414.1.fasta Affe30_trimmed_R1.fastq.gz Affe30_trimmed_R2.fastq.gz > Affe30_reads_aligned.sam Convert the SAM file to BAM and sort it: samtools view -Sb Affe30_reads_aligned.sam | samtools sort -o Affe30_reads_aligned_sorted.bam samtools index Affe30_reads_aligned_sorted.bam #20556 + 0 in total (QC-passed reads + QC-failed reads) #0 + 0 secondary #42 + 0 supplementary #0 + 0 duplicates #13645 + 0 mapped (66.38% : N/A) #20514 + 0 paired in sequencing #10257 + 0 read1 #10257 + 0 read2 #12582 + 0 properly paired (61.33% : N/A) #12648 + 0 with itself and mate mapped #955 + 0 singletons (4.66% : N/A) #0 + 0 with mate mapped to a different chr #0 + 0 with mate mapped to a different chr (mapQ>=5) h. Generate Coverage Profile for Contigs: The file Affe30_output_sorted.bam is already aligned against the reference genome. You can directly use it to visualize contig coverage. Ensure the BAM file is indexed: samtools index Affe30_contigs_aligned_sorted.bam i. Generate Coverage Tracks: Use BamCoverage to generate coverage files (in bigWig format) for both the reads and contigs. For reads coverage: bamCoverage -b Affe30_reads_aligned_sorted.bam -o Affe30_reads_coverage.bw For contigs coverage: bamCoverage -b Affe30_contigs_aligned_sorted.bam -o Affe30_contigs_coverage.bw j. calculate coverages samtools depth -a Affe30_reads_aligned_sorted.bam > coverage.txt awk '{sum+=$3} END {print "Average coverage = ",sum/NR}' coverage.txt #Average coverage = 26.6073 samtools depth -a Affe31_reads_aligned_sorted.bam > coverage.txt awk '{sum+=$3} END {print "Average coverage = ",sum/NR}' coverage.txt Average coverage = 22.0228 k. Visualize Alignments: Use tools like IGV (Integrative Genomics Viewer) or the following python scripts to visualize the alignment. import matplotlib.pyplot as plt import pandas as pd import pysam # File paths contig_bam_file = 'Affe30_contigs_aligned_sorted.bam' reads_bam_file = 'Affe30_reads_aligned_sorted.bam' # Function to calculate coverage from BAM file def calculate_coverage(bam_file): samfile = pysam.AlignmentFile(bam_file, "rb") coverage = [] # Iterate over each position in the BAM file to get coverage for pileupcolumn in samfile.pileup(): coverage.append(pileupcolumn.n) # Number of reads covering this position return coverage # Calculate read coverage read_coverage = calculate_coverage(reads_bam_file) # Create a DataFrame for read coverage read_positions = range(len(read_coverage)) # Position is the index read_data = pd.DataFrame({'position': read_positions, 'coverage': read_coverage}) # Plotting Read Coverage plt.figure(figsize=(12, 6)) plt.plot(read_data['position'], read_data['coverage'], color='blue') plt.title("Read Coverage Profile") plt.xlabel("Genomic Position") plt.ylabel("Coverage") plt.grid() plt.savefig("reads_coverage_profile.png") # Save the plot to a file plt.close() # Calculate contig coverage contig_coverage = calculate_coverage(contig_bam_file) # Create a DataFrame for contig coverage contig_positions = range(len(contig_coverage)) # Position is the index contig_data = pd.DataFrame({'position': contig_positions, 'coverage': contig_coverage}) # Plotting Contig Coverage plt.figure(figsize=(12, 6)) plt.plot(contig_data['position'], contig_data['coverage'], color='green') plt.title("Contigs Coverage Profile") plt.xlabel("Genomic Position") plt.ylabel("Coverage") plt.grid() plt.savefig("contigs_coverage_profile.png") # Save the plot to a file plt.close() print("Coverage plots saved successfully.") -
Reporting
cp ./Affe30_trimmed_vrap_out/spades/contigs.fasta Affe30_contigs.fasta cp ./Affe31_trimmed_vrap_out/spades/contigs.fasta Affe31_contigs.fasta cp ./Affe30_trimmed_vrap_out/blast/blastn.xlsx Affe30_contigs_annotation.xlsx cp ./Affe31_trimmed_vrap_out/blast/blastn.xlsx Affe31_contigs_annotation.xlsx 1. The contigs from the de novo assembly are incomplete and far from covering the full genome (see attached: Affe30_contigs.fasta, Affe31_contigs.fasta, Affe30_contigs_annotation.xlsx, Affe31_contigs_annotation.xlsx). 2. I've generated a plot showing both the reads and assembled contigs mapped to the closely related reference genome JX878414 (see attached: reads_and_contigs_on_JX878414.png). 3. One major issue I've noted is the low sequencing depth. After processing (trimming, adapter removal, and filtering out low-complexity reads), the remaining read count is very low. In comparison, the benchmark data from IDT shows they generated approximately 2-2.7 million paired reads per sample for a similar study. Our read count after quality control is significantly lower: Affe30 has only 10,257 paired-end reads, and Affe31 has 3,367 paired-end reads. Here's a summary of the read counts: Summary of Read Numbers * Raw read count: - Affe30: 88,183 x 2 paired-end reads - Affe31: 70,896 x 2 paired-end reads * After trimming (adapter sequences, bases with quality < Q30, and reads < 36 nt): - Affe30: 18,789 x 2 paired-end reads - Affe31: 11,282 x 2 paired-end reads * After removing low-complexity reads (e.g., GGGATTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTT): - Affe30: 10,257 x 2 paired-end reads - Affe31: 3,367 x 2 paired-end reads