fastqc -java.lang.OutOfMemoryError: Java heap space - Solved!

$fastqc filename.fastq
Started analysis of filename.fastq
Exception in thread "Thread-1" java.lang.OutOfMemoryError: Java heap space
    at uk.ac.babraham.FastQC.Modules.GCModel.GCModel.<init>(GCModel.java:69)
    at uk.ac.babraham.FastQC.Modules.PerSequenceGCContent.processSequence(PerSequenceGCContent.java:198)
    at uk.ac.babraham.FastQC.Analysis.AnalysisRunner.run(AnalysisRunner.java:88)
    at java.base/java.lang.Thread.run(Thread.java:844)


$ grep -n 'Xm' /home/prakki/sw/FastQC/fastqc
167:    unshift @java_args,"-Xmx${memory}m";
170:    unshift @java_args,'-Xmx250m';


$ vi /home/prakki/sw/FastQC/fastqc


Changed the 250m to 10g (depends on the RAM of your computer) in the line 170

$ fastqc filename.fastq
Started analysis of filename.fastq
Approx 5% complete for filename.fastq
Approx 10% complete for filename.fastq
Approx 15% complete for filename.fastq
Approx 20% complete for filename.fastq
Approx 25% complete for filename.fastq
Approx 30% complete for filename.fastq
Approx 35% complete for filename.fastq
Approx 40% complete for filename.fastq
Approx 45% complete for filename.fastq
Approx 50% complete for filename.fastq
Approx 55% complete for filename.fastq
Approx 60% complete for filename.fastq
Approx 65% complete for filename.fastq
Approx 70% complete for filename.fastq
Approx 75% complete for filename.fastq
Approx 80% complete for filename.fastq
Approx 85% complete for filename.fastq
Approx 90% complete for filename.fastq
Approx 95% complete for filename.fastq
Analysis complete for filename.fastq

Merged entries in the non-redundant Nucleotide (nt) database

I ran the following command for to extract a sequence from nt database:

$ blastdbcmd -db /mnt/Storage/nt/nt/nt -entry KJ413946.1 | grep '>'

>KJ413946.1 Escherichia coli strain ECS01 plasmid pNDM-ECS01, complete sequence >KP900017.1 Raoultella ornithinolytica strain YNKP001 plasmid pYNKP001-NDM, complete sequence

I noticed that this is a merged entry in the 'non-redundant' Nucleotide (nt) database, where two nucleotide sequences which have the same sequence are combined under a single FASTA entry.

When I used the second identifier in the merged entry, blastdbcmd command gave me the same result in return:

$ blastdbcmd -db /mnt/Storage/nt/nt/nt -entry KP900017.1 | grep '>'

>KJ413946.1 Escherichia coli strain ECS01 plasmid pNDM-ECS01, complete sequence >KP900017.1 Raoultella ornithinolytica strain YNKP001 plasmid pYNKP001-NDM, complete sequence

Curious to know, I checked if they are the same length. Yes they are of same length (41190).

$ blastdbcmd -db /mnt/Storage/nt/nt/nt -entry KP900017.1 | grep -v '>' | tr -d '\n' | wc 

      0       1   41190

$ blastdbcmd -db /mnt/Storage/nt/nt/nt -entry KJ413946.1 | grep -v '>' | tr -d '\n' | wc

      0       1   41190

Again curious to know,  if these two sequences are sharing an exact similarity match, I blasted both of them. And with no surprise, they share 100 % identity. See blast screenshots below:

So, how to sort this out? I need only one annotation instead of multiple for my analysis and do not want the second identifier to confuse my downstream analysis scripts.

By searching around, problem could solved by using  "-target_only"  option while running the command. I believe this situation might arise in nr database too but switching this on should help to achieve the purpose.

$ blastdbcmd -db /mnt/Storage/nt/nt/nt -entry KJ413946.1 -target_only | grep '>'

>KJ413946.1 Escherichia coli strain ECS01 plasmid pNDM-ECS01, complete sequence

Converting FastQ to FastA - comparison of one-liners and toolkits

Just out of curiosity, was checking time it takes to convert fastq to fasta by different tools and one-liners. I took the commands/one-liners from this biostar post.

The input file is gz file and contains 6 M reads of 301 bp.

Gubbins Phylogeny construction : Segmentation fault with some input files

 ***Note: There is a newer version gubbins (v2.3.4) which resolves the issue according to authors. Please find this post from the gubbins github page for more details.

This issue is already reported in gubbins github page for older versions, where it runs successfully for a set of fasta files and not for other set. Could not find much documentation to trace and resolve this error. This is the step which is throwing me the error:

gubbins -r -v seqnew.fa.gaps.vcf -a 100 -b 10000 -f test_gubbins/tmp/seqnew.fa -t seqnew.fa.iteration_1 -m 3 seqnew.fa.gaps.snp_sites.aln
Failed while running Gubbins. Please ensure you have enough free memory

After the parsnp alignment, i used the parsnp.fasta (seqnew.fasta here), to run the following steps. By breaking the alignment into shorter chunks I could successfully run gubbins. After the gubbins is run for all chunks, I combined them and plotted a phylogenetic tree. Commands I used are:

1) Formatting the alignment fasta file

$ sed -e '/^>/s/$/@/' -e 's/^>/#/' seqnew.fa | tr -d '\n'|tr "#" "\n"| tr "@" "\t" |sort -u -t ' ' -f -k 2,2 |sed '/^$/d'|sed -e 's/^/>/' -e 's/\t/\n/' >seqnew.oneline.fa ##fasta file in oneline    $ fasta_formatter -i seqnew.fa -o seqnew.oneline.fa -w 0 ##convert multiple fasta file to single line fasta file

2) Split the alignment file

$ perl splitAlignment.pl

-rw-rw-r-- 1 ttsh ttsh 7.7M Nov  9 11:00 seqnew.oneline.fa_0
-rw-rw-r-- 1 ttsh ttsh 7.7M Nov  9 11:00 seqnew.oneline.fa_1
-rw-rw-r-- 1 ttsh ttsh 7.7M Nov  9 11:00 seqnew.oneline.fa_2
-rw-rw-r-- 1 ttsh ttsh 7.7M Nov  9 11:00 seqnew.oneline.fa_3
-rw-rw-r-- 1 ttsh ttsh 7.7M Nov  9 11:00 seqnew.oneline.fa_4
-rw-rw-r-- 1 ttsh ttsh 7.7M Nov  9 11:00 seqnew.oneline.fa_5
-rw-rw-r-- 1 ttsh ttsh 2.8M Nov  9 11:00 seqnew.oneline.fa_6

splitAlignment.pl can be found at this github page.

3) Run Gubbins on chunk of sequences

$ for d in {0..6}; do time nohup run_gubbins.py --prefix postGubbins_ecloacae_$d --thread 12 --verbose --tree_builder fasttree seqnew.oneline.fa_$d >>gubbins.log 2>&1; done

4) Convert the alignment file into one line fasta format using fastx_toolkit

$ for d in {0..6}; do fasta_formatter -i postGubbins_ecloacae_$d.filtered_polymorphic_sites.fasta -o postGubbins_ecloacae_$d.filtered_polymorphic_sites.fasta_fastx -w 0; done

-rw-rw-r-- 1 ttsh ttsh 1.2M Nov  9 11:31 postGubbins_ecloacae_0.filtered_polymorphic_sites.fasta_fastx

-rw-rw-r-- 1 ttsh ttsh 1.3M Nov  9 11:31 postGubbins_ecloacae_1.filtered_polymorphic_sites.fasta_fastx
-rw-rw-r-- 1 ttsh ttsh 1.5M Nov  9 11:31 postGubbins_ecloacae_2.filtered_polymorphic_sites.fasta_fastx
-rw-rw-r-- 1 ttsh ttsh 1.3M Nov  9 11:31 postGubbins_ecloacae_3.filtered_polymorphic_sites.fasta_fastx
-rw-rw-r-- 1 ttsh ttsh 1.3M Nov  9 11:31 postGubbins_ecloacae_4.filtered_polymorphic_sites.fasta_fastx
-rw-rw-r-- 1 ttsh ttsh 1.2M Nov  9 11:31 postGubbins_ecloacae_5.filtered_polymorphic_sites.fasta_fastx
-rw-rw-r-- 1 ttsh ttsh 397K Nov  9 11:31 postGubbins_ecloacae_6.filtered_polymorphic_sites.fasta_fastx

5) Combine the alignment files

$ perl combinedAlignedfasta.pl 

combinedAlignedfasta.pl  can be found at this github page.

6) You can skip the 7th step by converting aligned fasta file to newick or .tre file using  phylogenetic tree inferring softwares like fasttree or RaxML.

$ FastTree -nt postGubbins_seqnew_ecloacae.filtered_polymorphic_sites.fasta > postGubbins_seqnew_ecloacae.filtered_polymorphic_sites.tre

6) Installed QIIME and from QIIME ran the make phylogeny script

$ qiime (enter into qiime environment)

$  qiime > make_phylogeny.py -i postGubbins_seqnew_ecloacae.filtered_polymorphic_sites.fasta -o postGubbins_seqnew_ecloacae.filtered_polymorphic_sites.tre

7) Used MEGA software to open the .tre file and plot the phylogenetic tree.

8) To change fasta file name to phyloID names, used a script taken from biostars.

$ while IFS= read -r aline; do filename=`echo "$aline" | cut -f1 -d "     " `; phyloID=`echo "$aline" | cut -f2 -d "  " `; sed -i "s/$filename/$phyloID/" postGubbins_seqnew_ecloacae.filtered_polymorphic_sites_Final_Parsed.nwk ; echo "$aline";  done < ~/Datta/7International_NDM/International_NDM_Proj_Filenames_PhyloIDs.txt 

$ perl script.pl PhyloIDs.txt postGubbins_seqnew_ecloacae.filtered_polymorphic_sites.tre >postGubbinsChunk.filtered_polymorphic_sitesv2.tre

Note: *" " in the above one liner of cut command are tabs

updated on June 29th, 2018

Code: hilite.me, lang: CSS, CSS: border:solid gray;border-width:.1em .1em .1em .8em;padding:.2em .6em;, Style: Perldoc

Converting NCBI annotation from GFF3 to BED12

boxed Language: Python, Style: perldoc

Recently, I needed to convert a gff3 file from NCBI to bed12 (bed file with 12 columns) format for UCSC browser purposes. The BED12 file is not available in UCSC since its latest genome is processed with the annotation. All, I needed to do was to convert gff3 to bed12. These are the steps involved:

Direct conversion (from GFF3 to BED12) using GFFtools 

1. (older version of script)

$ gff_to_bed.py ../ref_ASM185804v2_top_level.gff3 > ref_ASM185804v2_top_level_gff_to_bed_frm_Vipin.bed

$ head ref_ASM185804v2_top_level_gff_to_bed_frm_Vipin.bed
NC_031971.1 6018623 6018697 rna15273 . - 6018623 6018697 0 1 74, 0,
NW_017615921.1 193623 198461 rna67841 . - 193623 198461 0 6 633,85,89,39,44,68, 0,1075,1280,4477,4625,4770,
NW_017615921.1 193623 198461 rna67842 . - 193623 198461 0 7 633,85,89,354,39,44,68, 0,1075,1280,2684,4477,4625,4770,
NC_031968.1 33369464 33374965 rna7798 . - 33369464 33374965 0 2 1181,494, 0,5007,
NC_031980.1 33188366 33210198 rna43547 . + 33188366 33210198 0 4 916,183,9,4122, 0,12798,16099,17710,
NW_017615939.1 319272 321966 rna68222 . + 319272 321966 0 3 1350,253,749, 0,1526,1945,
NC_031976.1 23579551 23595921 rna31865 . + 23579551 23595921 0 26 115,113,87,49,163,120,258,731,455,528,155,83,94,389,156,236,250,218,133,141,155,97,167,350,147,443, 0,577,3241,4303,4441,4760,5592,6386,7302,7853,8486,10035,10364,10572,11048,11291,11621,11953,12301,13122,13352,13722,14853,15113,15669,15927,
NC_031976.1 23579551 23595921 rna31866 . + 23579551 23595921 0 26 115,113,87,49,163,120,258,731,455,528,155,83,94,389,156,236,250,218,133,141,155,97,167,350,147,395, 0,577,3241,4303,4441,4760,5592,6386,7302,7853,8486,10035,10364,10572,11048,11291,11621,11953,12301,13122,13352,13722,14853,15113,15669,15975,
NC_031976.1 23579551 23595921 rna31867 . + 23579551 23595921 0 27 115,113,87,49,163,120,258,731,138,227,528,155,83,94,389,156,236,250,218,133,141,155,97,167,350,147,443, 0,577,3241,4303,4441,4760,5592,6386,7302,7530,7853,8486,10035,10364,10572,11048,11291,11621,11953,12301,13122,13352,13722,14853,15113,15669,15927,
NC_031976.1 23579551 23595921 rna31868 . + 23579551 23595921 0 27 115,113,87,49,163,120,222,731,138,227,528,155,83,94,389,156,236,250,218,133,141,155,97,167,350,147,443, 0,577,3241,4303,4441,4760,5592,6386,7302,7530,7853,8486,10035,10364,10572,11048,11291,11621,11953,12301,13122,13352,13722,14853,15113,15669,15927,

$ wc -l  ref_ASM185804v2_top_level_gff_to_bed_frm_Vipin.bed
71103 ref_ASM185804v2_top_level_gff_to_bed_frm_Vipin.bed

2. (new version of script) - throws error

$ python gff_to_bed_latest.py ref_ASM185804v2_top_level.gff3 >out3.bed
Traceback (most recent call last):
  File "gff_to_bed_latest.py", line 112, in <module>
    __main__() 
  File "gff_to_bed_latest.py", line 109, in __main__
    writeBED(Transcriptdb)
  File "gff_to_bed_latest.py", line 66, in writeBED
    score = ent1['transcript_score'][idx] if ent1['transcript_score'].any() else score ## getting the transcript score 
ValueError: no field of name transcript_score

Raised an issue on the above error : https://github.com/vipints/GFFtools-GX/issues/9?_pjax=%23js-repo-pjax-container

Indirect conversion (from GFF3 to GTF to BED12)

1. Converting gff3 to genePred to bed12 - UCSC way

gff3ToGenePred - Convert a GFF3 file to a genePred file

$ ./gff3ToGenePred ../ref_ASM185804v2_top_level.gff3 ref_ASM185804v2_top_level.GenePred

genePredToBed - Convert from genePred to bed format.

$ ./genePredToBed ref_ASM185804v2_top_level.GenePred ref_ASM185804v2_top_level_genePredToBed.bed

This has generated me BED12 format successfully.

$ tail ref_ASM185804v2_top_level_genePredToBed.bed

NC_031965.1 434912 448252 rna12 0 - 436871 444514 0 5 3105,384,137,588,456, 0,9275,9991,10219,12884,
NC_031965.1 434912 448252 rna11 0 - 436871 444514 0 5 3105,384,137,502,456, 0,9275,9991,10219,12884,
NC_031965.1 339041 345377 rna10 0 - 341012 345350 0 2 94,5970, 0,366,
NC_031965.1 211588 212299 rna7 0 + 212021 212282 0 2 86,367, 0,344,
NC_031965.1 137999 139964 rna6 0 + 137999 139964 0 4 34,758,655,65, 0,137,1024,1900,

$ wc -l ref_ASM185804v2_top_level_genePredToBed.bed
63081 ref_ASM185804v2_top_level_genePredToBed.bed

2. Converting gff3 to gtf to bed12 

$ gffread ../ref_ASM185804v2_top_level.gff3 -T -o ref_ASM185804v2_top_level.gtf

$ ./cufflinks_gtf2bed.py ref_ASM185804v2_top_level.gtf >ref_ASM185804v2_top_level_cufflinks_gtf2bed.bed

$ head ref_ASM185804v2_top_level_cufflinks_gtf2bed.bed 
NC_013663.1 0 69 rna70978 100 + 0 69 255,0,0 1 69 0
NC_013663.1 69 1013 rna70979 100 + 69 1013 255,0,0 1 944 0
NC_013663.1 1013 1085 rna70980 100 + 1013 1085 255,0,0 1 72 0
NC_013663.1 1085 2784 rna70981 100 + 1085 2784 255,0,0 1 1699 0
NC_013663.1 2784 2858 rna70982 100 + 2784 2858 255,0,0 1 74 0
NC_013663.1 3836 3906 rna70983 100 + 3836 3906 255,0,0 1 70 0
NC_013663.1 3905 3976 rna70984 100 - 3905 3976 255,0,0 1 71 0
NC_013663.1 3976 4045 rna70985 100 + 3976 4045 255,0,0 1 69 0
NC_013663.1 5091 5163 rna70986 100 + 5091 5163 255,0,0 1 72 0
NC_013663.1 5164 5233 rna70987 100 - 5164 5233 255,0,0 1 69 0

$ wc -l ref_ASM185804v2_top_level_cufflinks_gtf2bed.bed
71912 ref_ASM185804v2_top_level_cufflinks_gtf2bed.bed

Now, all three methods generated three BED12 format files but the number of records are different. Let's check this by counting number of lines:

$ wc -l */*bed
   71103 direct/ref_ASM185804v2_top_level_gff_to_bed_frm_Vipin.bed
   71912 indirect/ref_ASM185804v2_top_level_cufflinks_gtf2bed.bed
   63081 indirect/ref_ASM185804v2_top_level_genePredToBed.bed

Number of records common to all three bed files (based on intervals):

$ cut -f2,3,4 */*bed | sort | uniq -c | sort -nk1,1 | fgrep ' 3 ' >common_to_threeBED_files_with_intervals.txt 
$ wc -l common_to_threeBED_files_with_intervals.txt
62070 common_to_threeBED_files_with_intervals.txt

Some records seem cleary missing using different conversion tools. I am not sure, if we can use these common 62070 features. But then, BED12 file could be generated from GFF3 file this way.

Observations:

1. gtf2bed threw following error. The gtf file was generated from gffread script
   
   

   $ gtf2bed < ref_ASM185804v2_top_level.gtf
   Error: Potentially missing gene or transcript ID from GTF attributes (malformed GTF at line [1]?)

   
   
2. Tried to eliminate the above error using the solution in biostar post from the author of script. Still without success.

$ awk '{ if ($0 ~ "transcript_id") print $0; else print $0" transcript_id \"\";"; }' ref_ASM185804v2_top_level.gtf | gtf2bed -
Error: Potentially missing gene or transcript ID from GTF attributes (malformed GTF at line [1]?)



3. Wanted to generate bed12 using this github script. Using the binary script gtfToGenePred from UCSC could not generate the required output

$ ./gtfToGenePred -genePredExt -geneNameAsName2 ref_ASM185804v2_top_level.gtf gene.tmp
ref_ASM185804v2_top_level.gtf doesn't appear to be a GTF file (GFF not supported by this program)

Reapr - Sampling fragment coverage - didn't get any coverage

 While running Reapr  on my genome using for the long insert libraries I received the following error.

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56

$ Reapr_1.0.18/reapr scaffolds.fasta LIMP_8kb_long_mapped.bam LIMP_8kb_Reapr_out perfect [REAPR pipeline] Running facheck [REAPR pipeline] Running preprocess [REAPR preprocess] sampling region scaffold1_size3169106:5458-6874. Already sampled 0 bases [REAPR preprocess] sampling region scaffold1_size3169106:6895-18716. Already sampled 1417 bases [REAPR preprocess] sampling region scaffold1_size3169106:30623-57811. Already sampled 13239 bases [REAPR preprocess] sampling region scaffold1_size3169106:57813-98842. Already sampled 40428 bases [REAPR preprocess] sampling region scaffold1_size3169106:98863-108371. Already sampled 81458 bases [REAPR preprocess] sampling region scaffold1_size3169106:108392-108786. Already sampled 90967 bases [REAPR preprocess] sampling region scaffold1_size3169106:108807-132036. Already sampled 91362 bases [REAPR preprocess] sampling region scaffold1_size3169106:153000-153621. Already sampled 114592 bases [REAPR preprocess] sampling region scaffold1_size3169106:153642-254312. Already sampled 115214 bases [REAPR preprocess] sampling region scaffold1_size3169106:254314-257386. Already sampled 215885 bases [REAPR preprocess] sampling region scaffold1_size3169106:257407-308464. Already sampled 218958 bases [REAPR preprocess] sampling region scaffold1_size3169106:309627-314663. Already sampled 270016 bases [REAPR preprocess] sampling region scaffold1_size3169106:314684-355810. Already sampled 275053 bases [REAPR preprocess] sampling region scaffold1_size3169106:355831-356388. Already sampled 316180 bases [REAPR preprocess] sampling region scaffold1_size3169106:356424-360553. Already sampled 316738 bases [REAPR preprocess] sampling region scaffold1_size3169106:362353-365556. Already sampled 320868 bases [REAPR preprocess] sampling region scaffold1_size3169106:365596-383718. Already sampled 324072 bases [REAPR preprocess] sampling region scaffold1_size3169106:393560-401076. Already sampled 342195 bases [REAPR preprocess] sampling region scaffold1_size3169106:401097-405358. Already sampled 349712 bases [REAPR preprocess] sampling region scaffold1_size3169106:405379-423388. Already sampled 353974 bases [REAPR preprocess] sampling region scaffold1_size3169106:423445-426298. Already sampled 371984 bases [REAPR preprocess] sampling region scaffold1_size3169106:426319-435419. Already sampled 374838 bases [REAPR preprocess] sampling region scaffold1_size3169106:435440-439265. Already sampled 383939 bases [REAPR preprocess] sampling region scaffold1_size3169106:439286-440463. Already sampled 387765 bases [REAPR preprocess] sampling region scaffold1_size3169106:440484-443034. Already sampled 388943 bases [REAPR preprocess] sampling region scaffold1_size3169106:443036-506052. Already sampled 391494 bases [REAPR preprocess] sampling region scaffold1_size3169106:506073-512313. Already sampled 454511 bases [REAPR preprocess] sampling region scaffold1_size3169106:512315-521665. Already sampled 460752 bases [REAPR preprocess] sampling region scaffold1_size3169106:526492-538394. Already sampled 470103 bases [REAPR preprocess] sampling region scaffold1_size3169106:538415-571424. Already sampled 482006 bases [REAPR preprocess] sampling region scaffold1_size3169106:571426-583083. Already sampled 515016 bases [REAPR preprocess] sampling region scaffold1_size3169106:583085-585679. Already sampled 526674 bases [REAPR preprocess] sampling region scaffold1_size3169106:585700-640937. Already sampled 529269 bases [REAPR preprocess] sampling region scaffold1_size3169106:640982-662781. Already sampled 584507 bases [REAPR preprocess] sampling region scaffold1_size3169106:662783-742983. Already sampled 606307 bases [REAPR preprocess] sampling region scaffold1_size3169106:743047-760078. Already sampled 686508 bases [REAPR preprocess] sampling region scaffold1_size3169106:765101-815949. Already sampled 703540 bases [REAPR preprocess] sampling region scaffold1_size3169106:815971-820741. Already sampled 754389 bases [REAPR preprocess] sampling region scaffold1_size3169106:820912-821293. Already sampled 759160 bases [REAPR preprocess] sampling region scaffold1_size3169106:821314-865486. Already sampled 759542 bases [REAPR preprocess] sampling region scaffold1_size3169106:865488-908573. Already sampled 803715 bases [REAPR preprocess] sampling region scaffold1_size3169106:909261-942195. Already sampled 846801 bases [REAPR preprocess] sampling region scaffold1_size3169106:942216-986026. Already sampled 879736 bases [REAPR preprocess] sampling region scaffold1_size3169106:986047-991833. Already sampled 923547 bases [REAPR preprocess] sampling region scaffold1_size3169106:991854-1003751. Already sampled 929334 bases [REAPR preprocess] sampling region scaffold1_size3169106:1003772-1027843. Already sampled 941232 bases [REAPR preprocess] sampling region scaffold1_size3169106:1029659-1034131. Already sampled 965304 bases [REAPR preprocess] sampling region scaffold1_size3169106:1034152-1034594. Already sampled 969777 bases [REAPR preprocess] sampling region scaffold1_size3169106:1034615-1057662. Already sampled 970220 bases [REAPR preprocess] sampling region scaffold1_size3169106:1057683-1063521. Already sampled 993268 bases [REAPR preprocess] sampling region scaffold1_size3169106:1063542-1064434. Already sampled 999107 bases [REAPR preprocess] Something went wrong sampling fragment coverage - didn't get any coverage. Mon Jul 18 11:42:50 SGT 2016

Solution that worked for me:

I changed the orientation of the mate pairs reads by reverse complementing before using REAPR since by default Reapr assumes that the orientation of the reads is FR. To perform RC on fastq files, I used fastx_reverse_complement. Check this page for more details.

Pilon Error - Inappropriate call if not paired read

I was running Pilon to improve my genome assembly using pilon-1.16.jar. 

time nohup java -jar /data/software/pilon-1.16.jar --genome genome.fasta --frags SIPE1_filtered_sorted.bam --frags SIPE2_filtered_sorted.bam --frags SIPE3_filtered_sorted.bam --frags SIPE4_filtered_sorted.bam --frags SIPE5_filtered_sorted.bam --frags SIPE6_filtered_sorted.bam --jumps LIMP1_filtered_sorted.bam --jumps LIMP4c_filtered_sorted.bam --output Pilon_Output --changes --threads 22 --verbose &

Pilon version 1.16 Mon Dec 7 09:45:45 2015 -0500

Genome: genome.fasta
Exception in thread "main" java.lang.reflect.InvocationTargetException
    at sun.reflect.NativeMethodAccessorImpl.invoke0(Native Method)
    at sun.reflect.NativeMethodAccessorImpl.invoke(NativeMethodAccessorImpl.java:62)
    at sun.reflect.DelegatingMethodAccessorImpl.invoke(DelegatingMethodAccessorImpl.java:43)
    at java.lang.reflect.Method.invoke(Method.java:498)
    at com.simontuffs.onejar.Boot.run(Boot.java:340)
    at com.simontuffs.onejar.Boot.main(Boot.java:166)
Caused by: java.lang.IllegalStateException: Inappropriate call if not paired read
    at htsjdk.samtools.SAMRecord.requireReadPaired(SAMRecord.java:648)
    at htsjdk.samtools.SAMRecord.getProperPairFlag(SAMRecord.java:656)
    at org.broadinstitute.pilon.BamFile$$anonfun$scan$1.apply(BamFile.scala:273)
    at org.broadinstitute.pilon.BamFile$$anonfun$scan$1.apply(BamFile.scala:268)
    at scala.collection.Iterator$class.foreach(Iterator.scala:750)
    at scala.collection.AbstractIterator.foreach(Iterator.scala:1202)
    at org.broadinstitute.pilon.BamFile.scan(BamFile.scala:268)
    at org.broadinstitute.pilon.GenomeFile$$anonfun$processRegions$4.apply(GenomeFile.scala:92)
    at org.broadinstitute.pilon.GenomeFile$$anonfun$processRegions$4.apply(GenomeFile.scala:92)
    at scala.collection.parallel.AugmentedIterableIterator$class.map2combiner(RemainsIterator.scala:115)
    at scala.collection.parallel.immutable.ParVector$ParVectorIterator.map2combiner(ParVector.scala:62)
    at scala.collection.parallel.ParIterableLike$Map.leaf(ParIterableLike.scala:1054)
    at scala.collection.parallel.Task$$anonfun$tryLeaf$1.apply$mcV$sp(Tasks.scala:49)
    at scala.collection.parallel.Task$$anonfun$tryLeaf$1.apply(Tasks.scala:48)
    at scala.collection.parallel.Task$$anonfun$tryLeaf$1.apply(Tasks.scala:48)
    at scala.collection.parallel.Task$class.tryLeaf(Tasks.scala:51)
    at scala.collection.parallel.ParIterableLike$Map.tryLeaf(ParIterableLike.scala:1051)
    at scala.collection.parallel.AdaptiveWorkStealingTasks$WrappedTask$class.internal(Tasks.scala:159)
    at scala.collection.parallel.AdaptiveWorkStealingForkJoinTasks$WrappedTask.internal(Tasks.scala:443)
    at scala.collection.parallel.AdaptiveWorkStealingTasks$WrappedTask$class.compute(Tasks.scala:149)
    at scala.collection.parallel.AdaptiveWorkStealingForkJoinTasks$WrappedTask.compute(Tasks.scala:443)
    at scala.concurrent.forkjoin.RecursiveAction.exec(RecursiveAction.java:160)
    at scala.concurrent.forkjoin.ForkJoinTask.doExec(ForkJoinTask.java:260)
    at scala.concurrent.forkjoin.ForkJoinPool$WorkQueue.runTask(ForkJoinPool.java:1339)
    at scala.concurrent.forkjoin.ForkJoinPool.runWorker(ForkJoinPool.java:1979)
    at scala.concurrent.forkjoin.ForkJoinWorkerThread.run(ForkJoinWorkerThread.java:107)

Elsewhere I read files with reads mapped in different directions is causing the problem. So, I tried Pilon1.16 with only short insert libraries direction, it worked for a bit of time but eventually failed!

Since, I wanted to run pilon using both paired end and mate-pair libraries, I tried downloading the newest version of Pilonv1.18 which ALSO DID NOT SOLVE the above error.

Solution that worked for me:

So after exchannging few mails with the author, I finally extracted the reads whose pairs also mapped elsewhere in the genome and generated a BAM file. This time alone Pilon worked and still running without error! Let us see..if it will last! :)

I used the following command from the samtools.

samtools view -b -f 0x2 LIMP1.bam | samtools sort - LIMP1_mappedPairs_sorted.bam