Skip to content

jielab/perhaps

Folders and files

NameName
Last commit message
Last commit date

Latest commit

 
 
 
 
 
 
 
 
 
 
 
 
 
 
 

Repository files navigation

PerHAPS (Paired-End Reads HAPlotyping for Sequencing)


#1. Download

  • Download perhaps_gui.exe and a few public tools (awk.exe, cut.exe, samtools.exe, sort.exe, uniq.exe ) into the "windows_tools" directory.

  • Download example data

id=NA20525
wget ftp.sra.ebi.ac.uk/vol1/run/ERR323/ERR3239807/$id.final.cram
samtools index $id.final.cram
samtools view -H XXX.cram | grep "SN:" | head -25 # check if the target BAM file XXX.cram has "chr" prefix

echo "1 159204012 159206500 ACKR1" > subset.bed
echo "19 44905781 44909393 APOE" >> subset.bed
sed -i 's/ /\t/g' subset.bed
samtools view -L subset.bed -O BAM -o $id.subset.bam $id.final.cram

#2. RUN

Linux version: Run "perhaps.sh" under "scripts" folder. Only the first 3 lines need to be updated.

Windows GUI version

  • Click perhaps_gui.exe to run the GUI version. The default value is pre-filled, and users only need to click the "submit" button.

  • The Windows version could also be called from the Windwos command terminal

python perhaps.py -i NA20525 -d .\test-data -s 1:159205564-159205704-159205737*

Figure 5

Figure 6

Figure 7

If users could not see the above images in browser, this is due to "DNS cache pollution". One short term fix for Windows users is to replace the "hosts" file (usually in "C:\Windows\System32\drivers\etc\hosts") with the "hosts" file posted on this site.


#3. Visualize

PerHAPS outputs a subset of the input BAM file that only contains paired short reads that are informative for the input haplotype.

Users could use IGV (http://www.igv.org/) to visualize the input BAM file and the output subset BAM file.

Figure 8


#4. Compare

make sure that VCF and BAM files are aligned to the same genome build

id=NA20525 # an example sample from G1K
  
# whatshap (https://whatshap.readthedocs.io/en/latest/)
  whatshap phase -o $id.whatshap.vcf --no-reference $id.vcf.gz $id.bam

# HapCUT2 (https://github.com/vibansal/HapCUT2)**	
  extractHAIRS --bam $id.bam --VCF $id.vcf --out $id.fragment
  HAPCUT2 --VCF $id.vcf --fragments $id.fragment --output $id.hap

# Smart-Phase (https://github.com/paulhager/smart-phase)
  java -jar smartPhase.jar -a $id.vcf.gz -p $id -g apoe.b38.bed -r $id.bam -m 60 -x -vcf -c 0.1 -o $id.tsv  


Contact & Cite

About

Paired End Reads based HAPlotyping for Sequencing (PMID: 33285565)

Resources

Stars

Watchers

Forks

Releases

No releases published

Packages

No packages published