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

One liner to convert a long fasta-file into many separate single fasta sequences

Taken from the Biostar post: https://www.biostars.org/p/105388/


while read line; do if [[ ${line:0:1} == '>' ]]; then outfile=${line#>}.fa; echo $line > $outfile; else echo $line >> $outfile; fi; done < combined_kpneumoniae.tfa

Coloring specific tips labels of phylogenetic tree in R

Some time ago, I wanted to change colors of specific branches and tips labels in a phylogenetic tree. I have R programming for changing colors. I have posted a screenshot of what I did.

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.

Replace a list of patterns with another list

$ while IFS= read -r aline; do filename=`echo "$aline" \| 

cut -f1 -d "        " `; phyloID=`echo "$aline" \| 

cut -f2 -d "  " `; 

sed -i "s/$filename/$phyloID/"  file_with_pattern_to_be_replaced.txt ; 

echo "$aline";  done < id_values.txt

## "id_values.txt" - contains pattern and replacement, separated by tab

## "file_with_pattern_to_be_replaced.txt" - replace patterns in this file

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)

Copying multiple files from a list to another directory

I have a list of filenames in a file "list.txt".

1. I want to find out the absolute location of these files
2. and then want to copy to another folder.

I did the following:

$ time for d in $(cat list.txt); do cp `find /media/My_Book/DATA/Raw/ -name "$d"` /media/F3F9/files/ ; done

Kernel driver not installed (rc=-1908) in Problem in Ubuntu 16.04 - Solved

1 2 3 4 5

$ sudo apt-get --reinstall install virtualbox virtualbox-dkms $ sudo apt-get install linux-headers-`uname -r` $ sudo apt-get install build-essential module-assistant $ sudo m-a prepare $ sudo modprobe vboxdrv

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.

Comparing bamcheck and samtools stats - finding inward and outward oriented read pairs from bam file

In addition to the biostar post, I am adding more details in this blog post.

I wanted to find out the inward and outward orientation pairs from a bam file. So, I checked a few tools. Two of them were bamcheck and samtools stats. The reason for doing so, if the bam shows that there are more number of outward oriented pairs, that means there are higher chances that the region has some problem associated with it.

But, I was stuck because both these tools gave different outputs. So, I myself wanted to find out which is one more reliable.

This subsequently led me run the following steps.

1) Split the bam into smaller bam files, (1 scaffold and its associated mapped reads)

$ bamtools split -in 500bp_insert_mapping_subset2_sorted.bam -reference

$ ls
500bp_insert_mapping_subset2_sorted.bam                   500bp_insert_mapping_subset2_sorted.REF__unitig_1900.bam
500bp_insert_mapping_subset2_sorted.REF__unitig_1023.bam  500bp_insert_mapping_subset2_sorted.REF__unitig_1933.bam
500bp_insert_mapping_subset2_sorted.REF__unitig_107.bam   500bp_insert_mapping_subset2_sorted.REF__unitig_1992.bam
500bp_insert_mapping_subset2_sorted.REF__unitig_1117.bam  500bp_insert_mapping_subset2_sorted.REF__unitig_2047.bam
500bp_insert_mapping_subset2_sorted.REF__unitig_1171.bam  500bp_insert_mapping_subset2_sorted.REF__unitig_2117.bam
500bp_insert_mapping_subset2_sorted.REF__unitig_118.bam   500bp_insert_mapping_subset2_sorted.REF__unitig_2142.bam
500bp_insert_mapping_subset2_sorted.REF__unitig_1226.bam  500bp_insert_mapping_subset2_sorted.REF__unitig_2151.bam
500bp_insert_mapping_subset2_sorted.REF__unitig_134.bam   500bp_insert_mapping_subset2_sorted.REF__unitig_2163.bam
500bp_insert_mapping_subset2_sorted.REF__unitig_1403.bam  500bp_insert_mapping_subset2_sorted.REF__unitig_2214.bam
500bp_insert_mapping_subset2_sorted.REF__unitig_1486.bam  500bp_insert_mapping_subset2_sorted.REF__unitig_228.bam
500bp_insert_mapping_subset2_sorted.REF__unitig_1527.bam  500bp_insert_mapping_subset2_sorted.REF__unitig_239.bam
500bp_insert_mapping_subset2_sorted.REF__unitig_1589.bam  500bp_insert_mapping_subset2_sorted.REF__unitig_312.bam
500bp_insert_mapping_subset2_sorted.REF__unitig_1630.bam  500bp_insert_mapping_subset2_sorted.REF__unitig_315.bam
500bp_insert_mapping_subset2_sorted.REF__unitig_1686.bam  500bp_insert_mapping_subset2_sorted.REF__unitig_363.bam
500bp_insert_mapping_subset2_sorted.REF__unitig_1812.bam  500bp_insert_mapping_subset2_sorted.REF__unitig_992.bam
500bp_insert_mapping_subset2_sorted.REF__unitig_1853.bam

2) Count the inward and oriented pairs, using both bamcheck 

$ for i in *bam; do echo "$i"; bamcheck $i | egrep '^SN.*inward|^SN.*outward'; done | paste - - - | column -t
500bp_insert_mapping_subset2_sorted.bam                   SN  inward  oriented  pairs:  95506  SN  outward  oriented  pairs:  0
500bp_insert_mapping_subset2_sorted.REF__unitig_1023.bam  SN  inward  oriented  pairs:  2801   SN  outward  oriented  pairs:  0
500bp_insert_mapping_subset2_sorted.REF__unitig_107.bam   SN  inward  oriented  pairs:  2861   SN  outward  oriented  pairs:  0
500bp_insert_mapping_subset2_sorted.REF__unitig_1117.bam  SN  inward  oriented  pairs:  2473   SN  outward  oriented  pairs:  0
500bp_insert_mapping_subset2_sorted.REF__unitig_1171.bam  SN  inward  oriented  pairs:  2844   SN  outward  oriented  pairs:  0
500bp_insert_mapping_subset2_sorted.REF__unitig_118.bam   SN  inward  oriented  pairs:  2620   SN  outward  oriented  pairs:  0
500bp_insert_mapping_subset2_sorted.REF__unitig_1226.bam  SN  inward  oriented  pairs:  3212   SN  outward  oriented  pairs:  0
500bp_insert_mapping_subset2_sorted.REF__unitig_134.bam   SN  inward  oriented  pairs:  3543   SN  outward  oriented  pairs:  0
500bp_insert_mapping_subset2_sorted.REF__unitig_1403.bam  SN  inward  oriented  pairs:  3018   SN  outward  oriented  pairs:  0
500bp_insert_mapping_subset2_sorted.REF__unitig_1486.bam  SN  inward  oriented  pairs:  3671   SN  outward  oriented  pairs:  0
500bp_insert_mapping_subset2_sorted.REF__unitig_1527.bam  SN  inward  oriented  pairs:  3530   SN  outward  oriented  pairs:  0
500bp_insert_mapping_subset2_sorted.REF__unitig_1589.bam  SN  inward  oriented  pairs:  3093   SN  outward  oriented  pairs:  0
500bp_insert_mapping_subset2_sorted.REF__unitig_1630.bam  SN  inward  oriented  pairs:  2951   SN  outward  oriented  pairs:  0
500bp_insert_mapping_subset2_sorted.REF__unitig_1686.bam  SN  inward  oriented  pairs:  1545   SN  outward  oriented  pairs:  0
500bp_insert_mapping_subset2_sorted.REF__unitig_1812.bam  SN  inward  oriented  pairs:  3242   SN  outward  oriented  pairs:  0
500bp_insert_mapping_subset2_sorted.REF__unitig_1853.bam  SN  inward  oriented  pairs:  2252   SN  outward  oriented  pairs:  0
500bp_insert_mapping_subset2_sorted.REF__unitig_1900.bam  SN  inward  oriented  pairs:  3876   SN  outward  oriented  pairs:  0
500bp_insert_mapping_subset2_sorted.REF__unitig_1933.bam  SN  inward  oriented  pairs:  3718   SN  outward  oriented  pairs:  0
500bp_insert_mapping_subset2_sorted.REF__unitig_1992.bam  SN  inward  oriented  pairs:  4269   SN  outward  oriented  pairs:  0
500bp_insert_mapping_subset2_sorted.REF__unitig_2047.bam  SN  inward  oriented  pairs:  2493   SN  outward  oriented  pairs:  0
500bp_insert_mapping_subset2_sorted.REF__unitig_2117.bam  SN  inward  oriented  pairs:  3140   SN  outward  oriented  pairs:  0
500bp_insert_mapping_subset2_sorted.REF__unitig_2142.bam  SN  inward  oriented  pairs:  4837   SN  outward  oriented  pairs:  0
500bp_insert_mapping_subset2_sorted.REF__unitig_2151.bam  SN  inward  oriented  pairs:  6686   SN  outward  oriented  pairs:  0
500bp_insert_mapping_subset2_sorted.REF__unitig_2163.bam  SN  inward  oriented  pairs:  2668   SN  outward  oriented  pairs:  0
500bp_insert_mapping_subset2_sorted.REF__unitig_2214.bam  SN  inward  oriented  pairs:  3145   SN  outward  oriented  pairs:  0
500bp_insert_mapping_subset2_sorted.REF__unitig_228.bam   SN  inward  oriented  pairs:  2622   SN  outward  oriented  pairs:  0
500bp_insert_mapping_subset2_sorted.REF__unitig_239.bam   SN  inward  oriented  pairs:  2729   SN  outward  oriented  pairs:  0
500bp_insert_mapping_subset2_sorted.REF__unitig_312.bam   SN  inward  oriented  pairs:  2880   SN  outward  oriented  pairs:  0
500bp_insert_mapping_subset2_sorted.REF__unitig_315.bam   SN  inward  oriented  pairs:  2668   SN  outward  oriented  pairs:  0
500bp_insert_mapping_subset2_sorted.REF__unitig_363.bam   SN  inward  oriented  pairs:  4305   SN  outward  oriented  pairs:  0
500bp_insert_mapping_subset2_sorted.REF__unitig_992.bam   SN  inward  oriented  pairs:  1814   SN  outward  oriented  pairs:  0

2) Count the inward and oriented pairs, using samtools stats

$ for i in *bam; do echo "$i"; ~/Documents/samtools-1.2/samtools stats $i | egrep '^SN.*inward|^SN.*outward'; done | paste - - - | column -t
500bp_insert_mapping_subset2_sorted.bam                   SN  inward  oriented  pairs:  96131  SN  outward  oriented  pairs:  295
500bp_insert_mapping_subset2_sorted.REF__unitig_1023.bam  SN  inward  oriented  pairs:  2807   SN  outward  oriented  pairs:  0
500bp_insert_mapping_subset2_sorted.REF__unitig_107.bam   SN  inward  oriented  pairs:  2861   SN  outward  oriented  pairs:  0
500bp_insert_mapping_subset2_sorted.REF__unitig_1117.bam  SN  inward  oriented  pairs:  2474   SN  outward  oriented  pairs:  0
500bp_insert_mapping_subset2_sorted.REF__unitig_1171.bam  SN  inward  oriented  pairs:  2845   SN  outward  oriented  pairs:  0
500bp_insert_mapping_subset2_sorted.REF__unitig_118.bam   SN  inward  oriented  pairs:  2622   SN  outward  oriented  pairs:  0
500bp_insert_mapping_subset2_sorted.REF__unitig_1226.bam  SN  inward  oriented  pairs:  3221   SN  outward  oriented  pairs:  0
500bp_insert_mapping_subset2_sorted.REF__unitig_134.bam   SN  inward  oriented  pairs:  3553   SN  outward  oriented  pairs:  0
500bp_insert_mapping_subset2_sorted.REF__unitig_1403.bam  SN  inward  oriented  pairs:  3020   SN  outward  oriented  pairs:  0
500bp_insert_mapping_subset2_sorted.REF__unitig_1486.bam  SN  inward  oriented  pairs:  3677   SN  outward  oriented  pairs:  0
500bp_insert_mapping_subset2_sorted.REF__unitig_1527.bam  SN  inward  oriented  pairs:  3545   SN  outward  oriented  pairs:  0
500bp_insert_mapping_subset2_sorted.REF__unitig_1589.bam  SN  inward  oriented  pairs:  3097   SN  outward  oriented  pairs:  0
500bp_insert_mapping_subset2_sorted.REF__unitig_1630.bam  SN  inward  oriented  pairs:  2961   SN  outward  oriented  pairs:  0
500bp_insert_mapping_subset2_sorted.REF__unitig_1686.bam  SN  inward  oriented  pairs:  1548   SN  outward  oriented  pairs:  0
500bp_insert_mapping_subset2_sorted.REF__unitig_1812.bam  SN  inward  oriented  pairs:  3251   SN  outward  oriented  pairs:  0
500bp_insert_mapping_subset2_sorted.REF__unitig_1853.bam  SN  inward  oriented  pairs:  2272   SN  outward  oriented  pairs:  24
500bp_insert_mapping_subset2_sorted.REF__unitig_1900.bam  SN  inward  oriented  pairs:  3885   SN  outward  oriented  pairs:  0
500bp_insert_mapping_subset2_sorted.REF__unitig_1933.bam  SN  inward  oriented  pairs:  3727   SN  outward  oriented  pairs:  0
500bp_insert_mapping_subset2_sorted.REF__unitig_1992.bam  SN  inward  oriented  pairs:  4279   SN  outward  oriented  pairs:  0
500bp_insert_mapping_subset2_sorted.REF__unitig_2047.bam  SN  inward  oriented  pairs:  2850   SN  outward  oriented  pairs:  227
500bp_insert_mapping_subset2_sorted.REF__unitig_2117.bam  SN  inward  oriented  pairs:  3176   SN  outward  oriented  pairs:  14
500bp_insert_mapping_subset2_sorted.REF__unitig_2142.bam  SN  inward  oriented  pairs:  4841   SN  outward  oriented  pairs:  0
500bp_insert_mapping_subset2_sorted.REF__unitig_2151.bam  SN  inward  oriented  pairs:  6735   SN  outward  oriented  pairs:  24
500bp_insert_mapping_subset2_sorted.REF__unitig_2163.bam  SN  inward  oriented  pairs:  2679   SN  outward  oriented  pairs:  0
500bp_insert_mapping_subset2_sorted.REF__unitig_2214.bam  SN  inward  oriented  pairs:  3150   SN  outward  oriented  pairs:  0
500bp_insert_mapping_subset2_sorted.REF__unitig_228.bam   SN  inward  oriented  pairs:  2625   SN  outward  oriented  pairs:  0
500bp_insert_mapping_subset2_sorted.REF__unitig_239.bam   SN  inward  oriented  pairs:  2732   SN  outward  oriented  pairs:  0
500bp_insert_mapping_subset2_sorted.REF__unitig_312.bam   SN  inward  oriented  pairs:  2885   SN  outward  oriented  pairs:  0
500bp_insert_mapping_subset2_sorted.REF__unitig_315.bam   SN  inward  oriented  pairs:  2671   SN  outward  oriented  pairs:  0
500bp_insert_mapping_subset2_sorted.REF__unitig_363.bam   SN  inward  oriented  pairs:  4312   SN  outward  oriented  pairs:  0
500bp_insert_mapping_subset2_sorted.REF__unitig_992.bam   SN  inward  oriented  pairs:  1830   SN  outward  oriented  pairs:  6

3) I selected this small BAM file: 500bp_insert_mapping_subset2_sorted.REF__unitig_2151.bam

Bamcheck stats

$ bamcheck 500bp_insert_mapping_subset2_sorted.REF__unitig_2151.bam | head -28
# This file was produced by bamcheck (2012-04-24)
# The command line was:  bamcheck 500bp_insert_mapping_subset2_sorted.REF__unitig_2151.bam
# Summary Numbers. Use `grep ^SN | cut -f 2-` to extract this part.
SN sequences: 14606
SN is paired: 1
SN is sorted: 1
SN 1st fragments: 7843
SN last fragments: 6763
SN reads mapped: 14606
SN reads unmapped: 0
SN reads unpaired: 1234
SN reads paired: 13372
SN reads duplicated: 0
SN reads MQ0: 9217
SN total length: 1437208
SN bases mapped: 1437208
SN bases mapped (cigar): 1436380
SN bases trimmed: 0
SN bases duplicated: 0
SN mismatches: 0
SN error rate: 0.000000e+00
SN average length: 98
SN maximum length: 101
SN average quality: 36.8
SN insert size average: 497.2
SN insert size standard deviation: 29.1
SN inward oriented pairs: 6686
SN outward oriented pairs: 0

Samtools stats

$ samtools stats 500bp_insert_mapping_subset2_sorted.REF__unitig_2151_outties_filtered.bam | head -37
SN raw total sequences: 14606
SN filtered sequences: 0
SN sequences: 14606
SN is sorted: 1
SN 1st fragments: 7843
SN last fragments: 6763
SN reads mapped: 14606
SN reads mapped and paired: 13526 # paired-end technology bit set + both mates mapped
SN reads unmapped: 0
SN reads properly paired: 13372 # proper-pair bit set
SN reads paired: 13526 # paired-end technology bit set
SN reads duplicated: 0 # PCR or optical duplicate bit set
SN reads MQ0: 9217 # mapped and MQ=0
SN reads QC failed: 0
SN non-primary alignments: 0
SN total length: 1437208 # ignores clipping
SN bases mapped: 1437208 # ignores clipping
SN bases mapped (cigar): 1436380 # more accurate
SN bases trimmed: 0
SN bases duplicated: 0
SN mismatches: 0 # from NM fields
SN error rate: 0.000000e+00 # mismatches / bases mapped (cigar)
SN average length: 98
SN maximum length: 101
SN average quality: 36.8
SN insert size average: 491.7
SN insert size standard deviation: 29.4
SN inward oriented pairs: 6735
SN outward oriented pairs: 24

Observe here, (6735+24)*2 != 13,526 but equals to 13,518 paired reads. (check from step 7)

4) First take out all the mapped read pairs (irrespective of direction) from the bam and convert into sorted bed file. 

$ fgrep -f <(bedtools bamtobed -i 500bp_insert_mapping_subset2_sorted.REF__unitig_2151.bam \|
cut -f4 | cut -d "_" -f1 \|
sort | uniq -c | sort -n \| 
fgrep '      2' | sed 's/      2 //g') <(bedtools bamtobed -i 500bp_insert_mapping_subset2_sorted.REF__unitig_2151.bam) \|
sort -k4,4 >500bp_insert_mapping_subset2_sorted.REF__unitig_2151.bed

5) Checking the numbers of reads paired

$ cat 500bp_insert_mapping_subset2_sorted.REF__unitig_2151.bed \| 
awk '{print $4,$2,$3,$6}' \| 
paste - - | wc -l

6763 (reads paired: 13526) - This number is same as the number generated by samtools stats see above SN reads paired: 13526.

But the bamcheck was giving a different result (SN reads paired: 13372)

6) Now, from the read pairs, I want to know how many are 

1. Inward oriented pairs (--> <--)
2. Outward oriented pairs (<-- -->)
3. Same direction pairs (-->--> or <--<--)

7) At first I tried running the following command, which gave me slightly higher number than the samtools stats given number

$ cat 500bp_insert_mapping_subset2_sorted.REF__unitig_2151.bed \|
awk '{print $4,$2,$3,$6}' | paste - - \|
awk '{if ($2 > $6) print $5,$6,$7,$8,$1,$2,$3,$4; else print $0;}' \|
sed 's/-/<@@@/g' | sed 's/+/###>/g' | awk '{print $1,$2,$3,$5,$6,$7,$4$8}' \|
column -t | awk '{print $NF}' | sort | uniq -c

    27  <--()-->
   6736 -->()<--

(27 + 6736)*2 =13,526 paired reads.

But, samtools stats gives (6735+24 ) * 2 = 13518 reads are paired reads. I tried understand deeply and scrutinizing the bed file I found that 3 reads pairs are overlapping at the same position for outward oriented pairs and 1 overlapping read pair from the inward oriented pairs.

From overlapping, I mean the coordinates of both reads are mapping around the same point.

HISEQ:120:H0BJ3ADXX:2:2104:7360:21274_2:N:0:ACAGTG/1   8700   8798   HISEQ:120:H0BJ3ADXX:2:2104:7360:21274_2:N:0:ACAGTG/2   8700   8798   <-- -->
HISEQ:120:H0BJ3ADXX:1:1112:2349:75045_2:N:0:ACAGTG/1   8702   8796   HISEQ:120:H0BJ3ADXX:1:1112:2349:75045_2:N:0:ACAGTG/2   8702   8794   <-- -->
HISEQ:120:H0BJ3ADXX:2:1213:8095:95602_2:N:0:ACAGTG/1   8717   8814   HISEQ:120:H0BJ3ADXX:2:1213:8095:95602_2:N:0:ACAGTG/2   8717   8814   <-- -->

These were problematic ones!

8) So, I made slight change in the above command

$ cat 500bp_insert_mapping_subset2_sorted.REF__unitig_2151.bed \|
awk '{print $4,$2,$3,$6}' | paste - - \|
awk '{if ($2 > $6) print $5,$6,$7,$8,$1,$2,$3,$4; else if ( $2 == $6) next; else print $0;}' \|
sed 's/-/<--/g' | sed 's/+/-->/g' \|
awk '{print $1,$2,$3,$5,$6,$7,$4"()"$8}' \|
column -t | awk '{print $NF}' \|
sort | uniq -c

     24 <--()-->
   6735 -->()<--

which is same as samtools stats command.

I have checked these reads in a visualization tool called tablet. It shows the directions of the pairs when the bam file is loaded into the tool.

**Note:

The final script I used is not a one-liner to find the orientation of paired reads for bigger bam files, because other conditions have to be taken in consideration too. Such as overlapping pairs with different starting coordinates (in my script I used only same coordinates of both read pairs to report).