PP2PG(Phasing pipeline based on

这个软件好像都没有什么人用哎~
但是我今天用上了，为什么会用上呢，还不是在KB的组装那个三倍体的基因组，为了区分测序回来的那么那么多的reads都是属于谁家的！！！！
对于这个软件，需要一些依赖，不过基本上这些依赖我们早就用过了，也就不存在依赖不依赖的了吧~

需要的依赖环境

浅浅列举一下吧：做到这一步，这些软件还没有的话，有点说不过去吧？

1.  MUMmer Version:4.0.0beta2 ([https://github.com/mummer4/mummer](https://github.com/mummer4/mummer)) 
2.  Bismark Version:v0.22.1 ([https://github.com/FelixKrueger/Bismark](https://github.com/FelixKrueger/Bismark)) 
3.  minimap2 Version:2.17 ([https://github.com/lh3/minimap2](https://github.com/lh3/minimap2)) 
4.  HISAT2 Version:2.1.0 ([http://ccb.jhu.edu/software/hisat2/manual.shtml](http://ccb.jhu.edu/software/hisat2/manual.shtml)) 
5.  pysam Version:0.15.2 ([https://pypi.org/project/pysam/](https://pypi.org/project/pysam/)) 
6.  SAMtools Version:1.6 ([http://samtools.sourceforge.net/](http://samtools.sourceforge.net/)) 

Installation:

 git clone https://github.com/jia-wu-feng/PP2PG.git
 cd PP2PG
 chmod +X tool/phasing.py 

这是github上的一张流程介绍

image.png

Fig.1 The phasing pipeline based on two parental genomes (PP2PG). It mainly consists of four steps: (i) Mapping, (ii) SNP Calling, (iii) Phasing and (iv) Classification.
也就是包括四个步骤：（i）映射，（ii）SNP calling，（iii）评分和（iv）分类。

(i)mapping(四选一)

1.Iso-Seq /PicBio读取的类型使用 minimap2 软件进行对齐：

minimap2 -ax splice -uf --secondary=no -C5 -O6,24 -B4 --MD ref query > results.sam      

2.利用hisat2软件对RNA-Seq配对双末端读段的类型进行比对：

hisat2 --dta -k 1 -x genome_index -1 query_1 -2 query_2 -S out.sam

3.BS-Seq读取的类型使用Bismark软件进行对齐：

bismark --bowtie2 -genome genome_index -1 query_1 -2 query_2

4.Ribo-seq读取的类型使用Tophat2（Bowtie2）软件进行对齐：

tophat2 -N 2 -I 50000 -G gff_file/gtf_file -o prefix genome_index query

(ii) SNP Calling

nucmer --mum --maxgap=500 --mincluster=100 --prefix=g1_g2 g1_genome.fa g1_genome.fa

delta-filter -1 -q -r g1_g2.delta > g1_g2.filter.delta

show-snps -C -H -I -T -r -l g1_g2.filter.delta > g1_g2.snp

(iii) Scoring and Classification:

phasing.py --pt=<phase_type> --g1=<genome1_alignment> --g2=<genome2_alignment> 
           --snp=<snpfile> --gop1=<genome1_output_prefix> --gop2=<genome2_output_prefix>

4.Output file

Two genome-based alignment results file, containing 2 folders. One is the phasing results of g1 genome alignment, which includes: g1-synteny reads, g2-synteny reads, unknown reads, and reads that only mapped to the g1 genome (g1-only reads); the other is the phasing results of g2 genome alignment, which includes: g1-synteny reads, g2-synteny reads, unknown reads, and reads that only mapped to the g2 genome (g2-only reads).
两个基于基因组的比对结果文件，包含 2 个文件夹。 一个是g1基因组比对的阶段结果，其中包括：g1-synteny读段，g2-synteny读段，未知读段和仅映射到g1基因组的读段（仅g1读段）;另一种是G2基因组比对的阶段结果，其中包括：G1-Synteny读段，G2-Synteny读段，未知读段和仅映射到G2基因组的读段（仅G2读段）。

g1_genome/g1reads.bam --- g1-synteny reads mapped to g1 genome
g1_genome/g2reads.bam --- g2-synteny reads mapped to g1 genome
g1_genome/unknowreads.bam --- unknown reads mapped to g1 genome
g1_genome/g1onlyreads.bam --- reads that only mapped to the g1 genome (g1-only reads)

g2_genome/g1reads.bam --- g1-synteny reads mapped to g2 genome
g2_genome/g2reads.bam --- g2-synteny reads mapped to g2 genome
g2_genome/unknowreads.bam --- unknown reads mapped to g2 genome
g2_genome/g2onlyreads.bam --- reads that only mapped to the g2 genome (g2-only reads)

g1_g2_Phasing_Report.txt --- Final Phasing Reads Report
g1_g2_Phasing_Report.txt ---最终阶段阅读报告

Example

==========================
Final Phasing Reads Report
==========================

Phasing Data Type:      Iso-Seq (PacBio Isoform Sequence).

Note: When calculating the Separation Rate, the reads which are unmapped in two parental genomes are discarded.

Type of reads:  Read Counts
MH63-synteny reads: 129
ZS97-synteny reads: 129
Unknown reads:  47
MH63-only reads:    2
ZS97-only reads:    3

Separation Rate:        83.23%

还是看看自己是怎么用的吧~~~

#  先划分参考的序列，
minimap2 -a -t 40 -x asm10 -m 50 --secondary=no A.contig.mmi m64189e_221213_192604.ccs.fasta >A.sam
minimap2 -a -t 40 -x asm10 -m 50 --secondary=no BB.contig.mmi m64189e_221213_192604.ccs.fasta >BB.sam
#sam文件实在是太占内存了，还是转为bam吧
samtools sort -@ 20 -O bam -o A.sorted.bam A.sam
samtools sort -@ 20 -O bam -o BB.sorted.bam BB.sam
#我就有比对质量高的reads的id
awk '$5 > 55 {print $1}' A.sorted.bam | sort -u > A.sorted.id
#把A和BB的sorted的id先取个交集
cat A.sorted.id BB.sorted.id | sort | uniq -c | awk '$1 == 2{print $2}' > A_BB.JJ.id
#再来取个差集
sort A.sorted.id A_BB.JJ.id A_BB.JJ.id > A.uniq.id
##all.reads.id - A_and_BB.uniq.id(取差集)
sort all.reads.id A_and_BB.uniq.id A_and_BB.uniq.id | uniq -u > unkown.reads.id
#提取unkown的reads序列
seqtk subseq ../ccs.fasta unkown.reads.id > unkown.reads.fa

上面的和这个软件没有什么关系，我就单纯想写，下面才是这个软件怎么用~~~

#step1.Mapping
minimap2 -ax splice -uf --secondary=no -C5 -O6,24 -B4 --MD -t 60 ../A.contig.fa unkown.reads.fa |samtools sort -@ 60 -O bam -o A.sorted.bam -
minimap2 -ax splice -uf --secondary=no -C5 -O6,24 -B4 --MD -t 60 ../BB.contig.fa unkown.reads.fa |samtools sort -@ 60 -O bam -o B.sorted.bam -

#Step2.SNP calling
nucmer --mum --maxgap=500 --mincluster=100 --prefix=A_BB -t 60 ../A.contig.fa ./BB.contig.fa
delta-filter -1 -q -r A_BB.delta > A_BB.filter.delta
show-snps -C -H -I -T -r -l A_BB.filter.delta > A_BB.snp

#step3.Scoring and Classification
/share/home/off/biosoft/PP2PG/tool/phasing.py --pt=Isoseq --g1=A.sorted.bam --g2=B.sorted.bam --snp=A_BB.snp --gop1=A.reads  --gop2=B.reads --st 60

我也在等结果！！！

参考资料

https://github.com/jia-wu-feng/PP2PG