Does VDJer work on BAMs from HISAT2?

[ooondrej]

Hi, I am trying to use VDJer to pull out BCR-seq data from sorted and indexed BAM files created by HISAT2 (alignment with hg38 ucsc). However, when I run it, it stops at PRE_PRE_GRAPH1 stage with the following error:

input_len: 7333, record_len: 151 Initial char in input invalid: J ERROR! INVALID INPUT: ===========================0TCATTCTTCTTTCAATCAGCAGGGACCGTGCACTCTCTTGGAGCCACCATAGAAAACAGAGGTGCATCCAGCACCACAGAAAACAGAGCCACCACAGAAAACAGAGGGTGACTAAFFFJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJFFJJJJF-FJJFJJJJJJJAFJJAJAJFAJAJJFJJFJJJJJJJJJJ<JJJJJJFJFJJJJF<<FAAF-FA0GCTGCAGCTGGGAGTGTGCAGAGACTGGAGGGGATGACAGTCACCCTCTGTTTTCTGTGGTGGCTCTGTTTTCTGTGGTGCTGGATGCACCTCTGTTTTCTATGGTGGCTCCAAF7F<--A)JJF-JJJJ<<<FFFFAFJJFAJJAFAJFAAF-FAAF<<FJJJJFJFJJJJJJ<JJJJJJJJJJFJJFJJAJAFJAJAJJFAJJJJJJJFJJF-FJJJJFFJJJJ0GTGACTGTCATCCCCTCCAGTCTCTGCACACTCCCAGCTGC ....

Does anyone know if VDJer accepts BAMs from aligners other than STAR? Could this error be because of the use of HISAT2 instead of STAR? Many thanks for any comments or help! Ondrej

[mozack]

Hi, It is possible that V'DJer can work on HISAT2 bam files, however we have not explicitly tested this.

The issue you are running into here is appears to be due to variable read lengths. As currently implemented, V'DJer requires all reads to be the same length.

[ooondrej]

Hi mozack, many thanks for this info! Can I just double check a few more related things, pls?

1. Does it mean that the reads in each pair have to be the same length, or all the reads in the fasq file have to be the same length? That would mean to apply a very strict trimming at the very beginning to achieve a uniform reads length.

2. Is the current default for --rl is 50? What is this parameter actually doing? Choosing only reads with a defined read length for the donwstream analysis from BAM? My fasq files contain reads about 140 bp long. Would it be reasonable for me to change --rl to higher number then?

Thanks a lot!

[ooondrej]

I have trimmed my sample fasq to a uniform reads length 100 across both paired files, aligned with HISAT2 and VDJer performed without any error message. However, my vdj_contigs.fa is empty and there is no sam file generated. This was bulk RNA-seq on sorted B cells.

command: vdjer --in test_sorted.bam --rl 100 --ins 175 --chain IGH --ref-dir /software/team152/vdjer-0.12/igh

Any thoughts what might be going wrong? Many thanks!

[mozack]

Can you confirm that the unmapped reads were output within the BAM file and also that hardclipping is not used for primary alignments?