Mapping Reads to Known Gene Paralogs with Long Read Technology

Question

I have some sequencing data from a captured region that is a known paralog edited. For now, I have been mapping the data using standard minimap2 flags for PacBio DNA sequencing:

minimap2 -ax map-pb GRCh37.fa fastqfile.fastq > samfile.sam

I captured a 10kb region and what I am seeing in IGV is puzzling. The region directly before and directly after the paralogues are covered at about 1000X, but the genes themselves are not covered at all:

The two genes in the middle are the paralogues. I left them off to protect data. This data is also only the sorted data, the duplicate marked data has lower coverage (600X).

Question 1

I was wondering if there is a different approach to mapping I could take that would output all alignments for a read, and then specifically only multi-mapped reads?

Question 2

Is this process normally accomplished by minimap2, or are other aligners more effective for paralogous mapping?

Question 3

The interesting thing is, while my average read length was 6kb, I am seeing the reads here are about 4kb. I have another region with very high coverage where the reads spanned both paralogues and mapping was successful.

Question: Why am I only seeing reads that are about 4kb long when the average read length for this dataset is 6kb?

Any help on this would be great! Thanks and I hope I have been clear on my question,

Dennis

Answer 1

Question 1

This command will get you all of supplementary alignments for the reads. This isn't exactly what you want though. You want all of the reads that have more than one mapping.

samtools view -f 2048 -h myfile_sorted.bam > supp_only.bam

This bash script returns a bam file that contains all alignments for reads that have multiple mappings. It works by only selecting reads that have supplementary alignments (2048 flag), or are not the primary alignment (256 flag), then using picard tools to filter the bam file. Using picard tools is the fastest way I have seen to filter a bam file given a list of reads.

#/bin/bash
BAMFILE=/path/to/my_file.bam
OUTFILE=/path/to/multiply_mapped.bam
PICARD_JAR=/path/to/picard.jar
TEMP=multiply_mapped.txt

for filt_int in 256 2048; \
do samtools view -f "$filtint" ${BAMFILE} | \
  grep -v '#' | cut -f1; done | sort | uniq > ${TEMP}; \
java -jar ${PICARD_JAR} FilterSamReads \
  I=${BAMFILE} \
  O=${OUTFILE} \
  READ_LIST_FILE=${TEMP} \
  FILTER=includeReadList
rm ${TEMP}

This command took 21 seconds to execute for my bam file with ~300,000 long reads mapped and several repetitive regions.

Question 2

Filtering the bam file should be pretty standard, as addressed above in my response to Q1. However, different long-read aligners are more sensitive to distantly-related paralogs. You might consider trying the blasr aligner since it is tuned for PacBio long reads, or try bwa-mem to compare your results with minimap2. Realistically though, most of these aligners will probably behave similarly to one another.

The fact that there are paralogs does not matter to the aligners. The aligner just sees these regions as two different places to which a single read can map just about equally well. When there are multiple similar regions to which a read can map, as in the case of paralogs, the aligner will report both alignments. I do not know of any long-read paralog-aware aligners.

As to why none of the reads mapped to the paralogs: Perhaps they have low mappability, as in they are repetitive? Also, maybe the capture didn't work properly? It looks like there is serious sequencing bias to the dataset given the uneven coverage.

Question 3

Are the reads that are mapping actually 4kb, or are the alignments only 4kb long? You may have many 6kb-long sequences that have been clipped at either the 5' or 3' end. If this is not a heavily-curated genome then the assembly may be incorrect around the paralogs, preventing proper mapping.

One suggestion is to take whole-genome long reads from the same organism and to map them. Then check if the mapping to this region is equally strange compared to your capture data. This will give you a better idea of what caused the uneven coverage and seeming lack of coverage around the paralogs.

Answer 2

I have found an alternative solution that has been excellent for this paralog mapping issue, so I thought I'd post what I had worked on:

Question 1+2

I was wondering if there is a different approach to mapping I could take that would output all alignments for a read, and then specifically only multi-mapped reads?

Using blasr, pacbio specified some options in their github docs for long, nearly identical reads:

blasr input.bam /path/to/reference.fa --hitpolicy all --bam --out alignments.bam --minMatch 8 --maxMatch 15 --nproc 4

The link to the docs as to why this is useful is here, although I did tweak the parameters to smaller seed sizes to create more anchors to the reference, which produced more secondary alignments. This is effective because this will output up to 10 secondary alignments per read, which I saw quite frequently for my region of interest.

After this, sorting was accomplished using Picard SortSam, with the following flags:

java -jar picard.jar SortSam INPUT=${SAMPLE} OUTPUT=./sorted/${SAMPLE}\
SORT_ORDER=coordinate VALIDATION_STRINGENCY=LENIENT\
CREATE_INDEX=true MAX_RECORDS_IN_RAM=100000

The blasr files have much more information than normal sams produced by aligners like minimap2 and bwa mem, so the MAX_RECORDS_IN_RAM=100000 as opposed to 500000 helps to ensure less frequent job failures (this was done on a cluster).

Next, I parsed out the region of interest:

samtools view -h input.bam "chrX:1000-2000" > extracted.bam # this is fake

Parse out mapq=0 reads:

awk '$5==0' extracted.bam > mapq0.txt

Get read names to grep for later:

awk '{print $1}' mapq0.txt > read_names_mapq0.txt

Search for all reads with readname from the original bam to get mapping positions:

samtools view infile.bam | grep -f read_names_mapq0.txt > all_reads_mapq0_mapq_high.txt

Keep only columns: 1) read_name, 2) Sam Flag, 3) Chr, 4) Pos, 5) MapQ:

awk '{print $1,$2,$3,$4,$5}' OFS='\t' all_reads_mapq0_mapq_high.txt > mapq_read_info.txt 

Then I wrote a simple python script to sort by read name, count the number of times a read shows up, and see all the positions it maps to. An example of the output:

        Read_Name       Flag    Chr     Pos     MAPQ
read1/10158955/ccs       272     19      3395674 0
read1/10158955/ccs       0       16      217682  254
read1/10158955/ccs       272     12      64990527        0
read1/10158955/ccs       256     16      222913  0
read1/10158955/ccs       272     19      21328060        0
read1/10158955/ccs       256     19      3224052 0
read1/10158955/ccs       272     17      34878621        0
read1/10158955/ccs       272     16      84204312        0
read1/10158955/ccs       256     1       244788680       0
read1/10158955/ccs       272     2       222951342       0
read2/10223868/ccs       16      16      217057  254
read2/10223868/ccs       256     20      3078724 0
read2/10223868/ccs       272     6       133521134       0
read2/10223868/ccs       256     7       152515722       0
read2/10223868/ccs       256     21      45258117        0
read2/10223868/ccs       256     1       249160685       0
read2/10223868/ccs       272     16      223631  0
read3/10551872/ccs       272     GL000220.1      144982  0
read3/10551872/ccs       256     18      72831954        0
...

read_counts output:

        Read_Name       Flag    Chr     Pos     MAPQ
read1/10158955/ccs       10      10      10      10
read2/10223868/ccs       7       7       7       7
read3/10551872/ccs       9       9       9       9

The numbers in each column correspond to the number of supplementary reads seen. A brief explanation of sam flags which I feel is relevant for unpaired reads:

Flag=0: read was mapped to + strand
Flag=16: read was mapped to - strand
Flag=256: read was mapped to + as secondary alignment
Flag=272: read was mapped to - as secondary alignment