0
$\begingroup$

I am working with a high performance cluster computer containing 112 threads. I am trying to align PacBio transcriptome reads against the genome to count the gene number. For pair end read i used the following workflow:

#convert gff to gtf
/home/software/cufflinks-2.2.1/gffread xxx.gff -T -o xxx.gtf
#build index
/home/software/hisat2-2.2.1/hisat2_extract_exons.py xxx.gtf > xxx.exon
/home/software/hisat2-2.2.1/hisat2_extract_splice_sites.py xxx.gtf > xxx.ss
/home/software/hisat2-2.2.1/hisat2-build -p 20 xxx.fa --ss xxx.ss --exon xxx.exon xxx_tran
#align reads
/home/software/hisat2-2.2.1/hisat2 -p 80 -x ./xxx_tran -1 R1.fq -2 R2.fq -S gonad.sam

#convert to bam and sort by read name
samtools view -bS gonad.sam -o gonad.bam

samtools sort -@ 80 gonad.bam -o gonad.sorted.bam

#count gene
htseq-count -f bam --strand=no gonad.sorted.bam xxx.all.gtf > gonad-counts.txt

But in case of Long read i am not sure what i need to do during the alignment stage. I have tried with minimap2 like as following:

#align reads
/home/software/minimap2-2.17_x64-linux/minimap2 -ax splice:hq -uf yyy.fa gill.fq > gill.sam

in following, i stucked in counting stage:

#count gene
htseq-count -f bam --strand=no gill.sorted.bam yyy.gtf > gill-counts.txt
error txt:
[E::idx_find_and_load] Could not retrieve index file for 'gill.sorted.bam'
100000 GFF lines processed.
200000 GFF lines processed.
300000 GFF lines processed.
400000 GFF lines processed.
438395 GFF lines processed.
[E::idx_find_and_load] Could not retrieve index file for 'gill.sorted.bam'
Error occured when processing input (record #79 in file gill.sorted.bam):
  'NoneType' object has no attribute 'encode'
  [Exception type: AttributeError, raised in _HTSeq.pyx:1379]

My question is that during counting it requires sorted bam and genome gtf file that was generated succesfully during previous staps, why it is getting failed even though the same pipeline is absolutely fine for pair end reads? Hisat2 can align long read? Is there any other way do that? As a novice i will be always grateful to you. Thanks.

Update: I tried generating an index:

# generate index
samtools index gill.sorted.bam

And the output was a .bai file like "gill.sorted.bam.bai". And I tried again to count as written above. But then I found the following error:

100000 GFF lines processed.
200000 GFF lines processed.
300000 GFF lines processed.
400000 GFF lines processed.
438395 GFF lines processed.
Error occured when processing input (record #79 in file gill.sorted.bam):
  'NoneType' object has no attribute 'encode'
  [Exception type: AttributeError, raised in _HTSeq.pyx:1379]

Is there any way?

Update 2:

Here is what the BAM file looks like:

SRR6418045.89565    2048    PYXC01000001.1  9893    60  5794H11M1I42M1I3M2I19M1D20M1I15M1D20M3D43M1I11M3D8M1I19M3I24M1I7M1I4M1I5M1D6M1D14M1I45M1D2M1D3M1D4M1I8M1I25M1D21M1D62M1D3M1D25M2D21M2I30M1I25M1D2M1D5M2I4M1D13M1D12M1D10M1I59M21564H    *   0   0   GAATACATAAACTGGATAGCGCACTGTTAATATTTACAGTGAAGAGAGTCTCTTCGAATTTTCCTCCTCAATCATCTTATTAGAAAGATGTTGAGTGGGCAGCAGCCAAATAAAATACTGTTCAGTGATCACTTCCAATGATGCTGTGCACTTACAGAAGCATGTATGGTACACAGTGTTATGCCTAATTTAAAAAACTAATGTAAATTATTTAAAAGTTTTTTTTTTTTTTATGAACTTTTTAGCCACATTGAATGGCTTCAGAATATCTACAGCCTAATGTATTTTCTTAGTAGCCTAATGAAGAAAGAATATTACCTTAAGTTACCATAGATTGACGCCAAAACTCGTCAAATCTCAAACCCCTAACATTTGTTCAACTGTTGTATACATGCTCATACACAGTGTTTCTTGTGTTCCTTAACCCACCGGGTAAAAAGTGGAGGTCCATGTCTGATACGTAGTTATCAATGTACTGTTGCCCATTTAGAGAGCACTTTGGTTTTCCCTCACCAGATAAGGTTTTTTCTCACTTTCAGTGTAAACCTGAGACTTGTTTGCACTTAATATTGGGGATTGTGTAACTCTTATTCAAGGAATTGCCTGTTTTCCAATTTGGTCTGATTACACTAACCTTTGTGTAAAATTATGGTTCAAATAAACTGTATCACA    -)/////,+$/"/..+//-/-+/.-/-.*,./.,)//'//.+/./..-/./../"/*/.),#)*.)-+)*-(///&.//.*'/,*////*-,./+.++#",.///,/,*,/-,)/,/)/.-)/.///./-..-*++)-(/*,../..)+//.////+*...////./...(././,/.,*.-'*/)(.-*)&(#'&,(#'....').,,$*')+#&,,)%%%$&,)*)*(*,..,'/-.+*)*,)+',-).-)*//,-",-+./)///'-.//,-+,.*//.-/%,**)-..../.*/,*'/.*//,*/.,(&.,+-)-'-*/*-/.-/-&+$&,'.'.&,.*,)/-."/,.,*/-//.,)/*,%.+*)/+.,--,./*--/.-.//'///)/--/$&//-//.//.-)/),,./-)+&../*,*,..-/+,)/,,,-*,/-,/.(.(///./(////./&-./-/+//*...////,&%+$&%#%$$--.././/.-+',-,#(.&%(-/,*/.///*..--++*)/-.*..*)*(,$*/--*.+./'//...-..)-.&.-)&))+.-,&#'-'--..'/'-//',//-*.)/+-*....%../.'(/*-$./,(.(////.#/..+/+%/'-'%...//*$+#--+/&-,-/,*-.-*/-.&../-,%/    NM:i:55 ms:i:343    AS:i:343    nn:i:0  ts:A:+  tp:A:P  cm:i:77 s1:i:460    s2:i:0  de:f:0.0656 SA:Z:PYXC01000001.1,9893,+,4147S672M2I23209S,60,37;PYXC01000001.1,9893,+,847S669M26I26488S,60,42;PYXC01000001.1,9893,+,26975S673M3I379S,60,44;PYXC01000001.1,9910,-,10980S654M9I16387S,60,41;PYXC01000001.1,9900,-,20736S661M9I6624S,60,43;PYXC01000001.1,9893,+,7409S673M15I19933S,60,49;PYXC01000001.1,9893,-,27284S667M12I67S,60,52;PYXC01000001.1,9893,-,19099S673M3I8255S,60,48;PYXC01000001.1,9893,-,7656S673M1652I18049S,4,1688;PYXC01000001.1,9908,+,2517S649M17I24847S,60,51;PYXC01000001.1,9893,+,23715S672M1D3643S,60,51;PYXC01000001.1,9893,-,25633S669M22I1706S,60,53;PYXC01000001.1,9893,-,17479S668M1I9882S,60,53;PYXC01000001.1,9893,-,14214S673M25I13118S,22,55;PYXC01000001.1,9899,+,13933S663M13I13421S,60,55;PYXC01000001.1,9900,-,22352S658M21I4999S,60,55;PYXC01000001.1,9901,+,15556S665M7I11802S,60,60;PYXC01000001.1,9893,+,12284S668M9I15069S,60,58;PYXC01000001.1,9909,-,24008S657M16I3349S,60,57;PYXC01000001.1,9902,+,25358S661M10I2001S,60,60;PYXC01000001.1,9893,+,10659S671M2D16700S,60,60;PYXC01000001.1,9893,+,9039S668M10I18313S,60,51;PYXC01000001.1,9893,-,2790S673M21I24546S,60,65;PYXC01000001.1,9899,+,22114S669M8I5239S,60,59;PYXC01000001.1,9893,-,6025S673M17I21315S,60,64;PYXC01000001.1,9908,-,12604S645M13I14768S,60,64;PYXC01000001.1,9893,-,4418S673M3I22936S,60,71;PYXC01000001.1,9893,-,1162S673M15I26180S,60,70;PYXC01000001.1,9898,+,17178S641M1I10210S,60,63;PYXC01000001.1,9893,+,20477S673M13I6867S,60,67;PYXC01000001.1,10390,-,64S160M2D27806S,60,8;   rl:i:2011
SRR6418045.89565    2048    PYXC01000001.1  9893    60  12284H12M1D27M1D9M1I15M1I26M1D19M2I16M2I2M1I17M1I32M1I31M2D8M3I25M1I6M1I19M2I12M1D32M2D15M1I6M2I7M1D4M2D17M1I10M1I32M1I36M1D4M1D2M1D20M1I59M1I11M1D10M1D10M1I24M1I9M2D2M1D28M1I24M1I11M15069H   *   0   0   GAATACATAAATGATAGCGCACTGTTAATATTTACAGTGAGAGAGTCTACTTGAATTCCTCCTCAAATCATCTTATTTAGAAAGATGTTGATGGGAGCAACCAAATAAAAATCTACTGTTCAGTGATCCAACTTTCACAAAATGATGCTGTTGCACTTACAGAAGCATGTATGGTACACAGTGTTATGACTATTTTAATAAAAAACAATGTAAATTTTAAAAGGTTTTTTTTTTTTTATGAACTTTTTAGCAACATTGAATGGTTCAGAAATATCTTAGCAAGCCTAATGTATTCTTAGTAGCCTAATGAAGAAAGAATATTACTAAGTTACCAGTAGGAATCTGAACGCCAAAATGTCATCTCAAACCCCTAACATTTTGGTTCAATCTGTTGTATACATGCATCATACACAGTGTTTCTTTGTGTTCCTTTAACCACCGGGTAAAAAGTGGAGGTCATGCTATTACTGTAGTTATCAATGTAACTGTTGCCCACTTTTAGAGAGCACTTTGGTTTTCTCACCAGATAAGGTTTTTTCTCACTTTTTATGTAAACTGAGACTTGTTGCACTGTACTATATTGAATATGTGTAACTCTTATTTTCAAGGATGCCTGTTTTCCATTTGGTCTATTAACACTAACACTTTGTGTAAACTTATGGTTCAAAATAAACTGT   -*//..//-*./%/,/./,,///.,/+//*/,////..,*..+/..-,(*.,.%&".*./%-..((//$///...)+-&.+'/////.+*$(),,../$.'.,)//-++$.."//,-.,/////'//,/$...*//./-()./+//,///,...%/-//,*.*/.//&/.//...///./..(-.///#/.+$(#.'..((%*(.)(-+-..%.*,((,,(..(-&+*)*******..+/$-..,++..--$-..'-.)//..-)/.,,*...//./.#.&/)*/)*(/./*,.&-(+////*/.).).+//-*//*)/*.*..-(//-/&,/-//..)*//.+.$/.',-+,(+(&//..//.,).%')./*///.--)--*.'(/"//...//.././/-)'(/,'..)%/./..&..$*,+-.-/*/++-%.+--+.),-,#,&(//.-,).-/./-.'./..///....-(*-/)/.#-$.//.%$*(*..())+/--/.../..(,+,..*).-/#.*'.+/-(/...&'+)/--.+'-))&$,'//*'/.-/(/.).,/,-,*.,/.*.".)..)./*/*/+/.'/%/&-),/..,)).(/--/&,*/..,*)-#+/.,-,./-**./%.//+/*)().(+-,,-.%'//../--.--/)*#..,*///.-   NM:i:58 ms:i:325    AS:i:325    nn:i:0  ts:A:+  tp:A:P  cm:i:65 s1:i:397    s
Improve this question
$\endgroup$
0

1 Answer 1

Reset to default
1
$\begingroup$

This error:

[E::idx_find_and_load] Could not retrieve index file for 'gill.sorted.bam'

Means that you have not generated an index for a BAM file. This needs to be done for most downstream processing of BAM files:

# generate index
samtools index gill.sorted.bam
# count genes
htseq-count -f bam --strand=no gill.sorted.bam yyy.gtf > gill-counts.txt

Update (after second error):

I've got no idea what's going wrong; this is not an error that I would expect to see, and doesn't look related to PacBio input. I wonder if the error is related to the number of threads used in the sort command, given that it happens at record 79, and you've specified 80 threads. The next things that I would do would be to look at the input files around where the error is to try to work out what's going wrong:

samtools view gill.sorted.bam | head -n 85 | tail -n 10
Improve this answer
$\endgroup$
2
  • $\begingroup$ Yes, i have used 80 threads. Do i need to reduce the thread number? $\endgroup$
    – Kishor Kumar Sarker
    Commented Aug 9, 2021 at 9:39
  • $\begingroup$ That depends on what machine you're running on. Attempting to run 80 threads on a machine that only has 4, for example, will slow down computation and may result in unexpected behaviour. $\endgroup$
    – gringer ♦
    Commented Aug 9, 2021 at 23:44

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.