Extracting strand specific reads from MinION cDNA-PCR protocol

Question

I recently performed my first MinION run, and now I'm trying to analyze the data. Being pretty new to the bioinformatics field, I was hoping some of you could help me out.

As a bit of background, I'm trying to define transcript isoforms, and down the line, I'm hoping to quantify isoform abundances. As my first attempt, I used a 1D kit that first performs RT starting at the polyA tail and then template switching to enrich for full length transcripts. To initiate RT, a poly-dT primer is added and after template switching a second primer is added. Using those added primer sequences, I then performed PCR to amplify the full-length cDNAs. This seemed to work well, so I went ahead and sequenced the library.

Now, on to my problem. I have aligned my fastq file to my genome with both bwa-mem and minimap2. In both cases, I get alignment of reads to the genome where I expect transcription to be occurring, but half of the reads are in the sense direction and half are in the antisense direction. This is presumably due to the fact that the cDNA is double stranded and either strand can be recruited to pore to be sequenced. However, because I'm interested in mapping transcript isoforms which could include anti-sense transcripts, this isn't exactly what I hoped for.

I'd like to figure out a way to extract only the reads that correspond to the actual mRNA transcript's sequence and not the reverse complement. Does anyone have ideas as far as programs that would help with this?
I'm thinking that I should be able to extract out only the molecules that have the polyA or have the strand switching primer sequence or something like that, but I'm not sure how to properly implement this.

Answer 1

This has been set up as a bioinformatics protocol on protocols.io.

I've been using LAST to identify the adapter orientation relative to the genome, and then using that information to split BAM files up and recombine them to create two strand-specific files that are displayable in a genome browser.

I have written my own script to process LAST results into a CSV format, which makes it easier to do line-by-line data filtering.

Note: these scripts have been slightly modified from scripts that I have run. Consider them demonstrative: pay attention to the words rather than the script.

Read correction

I prefer starting off my data analysis with a read correction with Canu (ideally v1.7 now that it's out, because that attempts correction of all reads, but here I use canu v1.6). I use minimap as the mapper to speed this up. The genomeSize parameter should be approximately a tenth to a fortieth of the number of bases in your dataset to make sure that no sequences are excluded (bigger is better, as long as Canu doesn't freak out about memory consumption):

~/install/canu/canu-1.6/Linux-amd64/bin/canu overlapper=minimap \
  genomeSize=100M minReadLength=100 minOverlapLength=30 -correct \
  -p 4T1_BC06 -d 4T1_BC06 -nanopore-raw \
  workspace/pass/barcode06/fastq_runid_*.fastq

This creates a file 4T1_BC06/4T1_BC06.correctedReads.fasta.gz.

Chimeric read filtering

The next step I carry out is a basic read-level QC to exclude chimeric reads. Porechop can be used for this, although that removes adapters by default, which is not particularly useful in this case.

I use LAST to search for adapter sequences within the corrected reads, pass it through my conversion script, and extract out duplicated mappings (i.e. where the same read/adapter pair appears more than once in the mapping results):

lastal -P 10 ONT_barcodes_adapters.fa <(zcat 4T1_BC06/4T1_BC06.correctedReads.fasta.gz) | \
  ~/scripts/maf_bcsplit.pl | awk -F',' '{print $1,$2}' | sort | \
  uniq -d | awk '{print $1}' | uniq > reads_with_duplicated_adapters.txt

I use this file to filter out chimeric reads from the corrected dataset using another fastq filtering script I've created:

~/scripts/fastx-fetch.pl -v -i reads_with_duplicated_adapters.txt 4T1_BC06/4T1_BC06.correctedReads.fasta.gz | \
  gzip > 4T1_BC06.correctedReads.uniqueOnly.fasta.gz

Adapter filtering

Along roughly similar lines to the chimeric read filtering, I then look for the strand switch and VNT adapters in the sequences. For a forward-orientated query, I expect the strand switch primer to be in the forward direction, and the VNP primer to be in the reverse direction:

lastal -P 10 ONT_barcodes_adapters.fa <(zcat 4T1_BC06.correctedReads.uniqueOnly.fasta.gz) | \
  ~/scripts/maf_bcsplit.pl | grep -e 'ONT_SSP,+' -e 'ONT_VNP,-' | \
  awk -F',' '{print $1}' | sort | uniq > fwdQry_seqs_BC06.txt
lastal -P 10 barcodes_primerSeqs.fa <(zcat 4T1_BC06.correctedReads.uniqueOnly.fasta.gz) | \
  ~/scripts/maf_bcsplit.pl | grep -e 'ONT_SSP,-' -e 'ONT_VNP,+' | \
  awk -F',' '{print $1}' | sort | uniq > revQry_seqs_BC06.txt

And then more filtering to split the reads up into forward and reverse [genome-relative] subsets:

~/scripts/fastx-fetch.pl -i fwdQry_seqs_BC06.txt 4T1_BC06.correctedReads.uniqueOnly.fasta.gz | \
  gzip > fwd_4T1_BC06.correctedReads.uniqueOnly.fasta.gz
~/scripts/fastx-fetch.pl -i revQry_seqs_BC06.txt 4T1_BC06.correctedReads.uniqueOnly.fasta.gz | \
  gzip > rev_4T1_BC06.correctedReads.uniqueOnly.fasta.gz

Mapping

Now that the reads have been oriented, the mapping can be done:

~/install/minimap2/minimap2 -t 10 -a mmus_ucsc_all_cdna.idx -x splice \
  fwd_4T1_BC06.correctedReads.uniqueOnly.fasta.gz \
  > fwd_4T1_BC06.CU_vs_mmus.bam
~/install/minimap2/minimap2 -t 10 -a mmus_ucsc_all_cdna.idx -x splice \
  rev_4T1_BC06.correctedReads.uniqueOnly.fasta.gz \
  > rev_4T1_BC06.CU_vs_mmus.bam

Splitting and recombining

The mapped files are then split based on their mapping direction using samtools view flag filtering to exclude or include reverse-mapped reads:

samtools view -F 0x10 -b fwd_4T1_BC06.CU_vs_mmus.bam > fwd_fwd_4T1_BC06.CU_vs_mmus.bam
samtools view -f 0x10 -b fwd_4T1_BC06.CU_vs_mmus.bam > rev_fwd_4T1_BC06.CU_vs_mmus.bam
samtools view -F 0x10 -b rev_4T1_BC06.CU_vs_mmus.bam > fwd_rev_4T1_BC06.CU_vs_mmus.bam
samtools view -f 0x10 -b rev_4T1_BC06.CU_vs_mmus.bam > rev_rev_4T1_BC06.CU_vs_mmus.bam

And finally the files are recombined to identify forward and reverse-encoded transcripts via samtools merge. I'm calling the "pos" strand the one that is encoded in the same direction as the genome, and the "neg" strand the reverse-complement direction:

samtools merge pos_4T1_BC06.CU_vs_mmus.bam \
  rev_fwd_4T1_BC06.CU_vs_mmus.bam fwd_rev_4T1_BC06.CU_vs_mmus.bam
samtools merge neg_4T1_BC06.CU_vs_mmus.bam \
  fwd_fwd_4T1_BC06.CU_vs_mmus.bam rev_rev_4T1_BC06.CU_vs_mmus.bam

Alternatively, it's possible to set the BAM flags for read 1 and read 2 for the reverse and forward-encoded adapters respectively, which allows everything to then be combined into a single BAM file (and treated in the same way as a strand-specific Illumina BAM file).

If this has worked properly, then mapping human or mouse to the mitochondrial genome should show most expression appearing on the positive strand, with a small scattering of negative-strand expression, a bit like this:

For a bit of a background on how/why this works, see this answer.