3.2.1.1 Convert genotype coding style

To take full advantage of matrix manipulation in R language, A/T/C/G coded genotypes will be first converted to 0/1 coded format with the base2num() function by providing reference and alternative alleles.

### covert A/T/C/G to 0/1
hetDa <- gmt[,1:4]
ref <- hetDa$ref
alt <- hetDa$alt

gmtDa <- gmt[,-(1:4)]
gmtDa <- base2num(gmt = gmtDa, ref = ref, alt = alt)
head(gmtDa)

##         gmt1 gmt2 gmt3 gmt4 gmt5
## 960055    NA   NA    1   NA   NA
## 1078815    0   NA   NA    0   NA
## 1079001    0   NA   NA    0    0
## 1079015    0   NA   NA    0    0
## 1079040   NA    1    1    0   NA
## 1079112   NA   NA    1    0   NA

3.2.1.2 Filter out genotyping error

hetSNPs with potential genotyping errors can be detected and filtered out by a Hidden MarKov Model (HMM) in the hapiErrorFitler() function. The transition and emission probabilities can be defined by users according to genotyping error rate and recombination frequency.

### define HMM probabilities
hmm = initHMM(States=c("S","D"), Symbols=c("s","d"), 
            transProbs=matrix(c(0.99999,0.00001,0.00001,0.99999),2),
            emissionProbs=matrix(c(0.99,0.01,0.01,0.99),2), 
            startProbs = c(0.5,0.5))
hmm

## $States
## [1] "S" "D"
## 
## $Symbols
## [1] "s" "d"
## 
## $startProbs
##   S   D 
## 0.5 0.5 
## 
## $transProbs
##     to
## from       S       D
##    S 0.99999 0.00001
##    D 0.00001 0.99999
## 
## $emissionProbs
##       symbols
## states    s    d
##      S 0.99 0.01
##      D 0.01 0.99

### filter out genotyping errors
gmtDa <- hapiFilterError(gmt = gmtDa, hmm = hmm)

## 35 hetSNPs with potential genotyping errors are filtered out !

3.2.1.3 Select high-quality hetSNPs to form a framework

For low coverage sequencing, hetSNPs that are genotyped in at least 3 gametes can be selected to form a ‘precursor’ framework for the draft haplotype construction. For dataset with high genotyping rate, hetSNPs that are observed in more gametes may be selected.

### select a subset of high-quality markers
gmtFrame <- hapiFrameSelection(gmt = gmtDa, n = 3) ###

## Number of hetSNPs in the framework: 21074

3.2.1.4 Imputation of missing data in the framework

Missing genotypes in each gamete are iteratively imputed by observed genotypes in other gametes to faciliate the draft haplotype inference. In the hapiImpute() function, nSPT specifies the minumum number of supporting gametes to impute a hetSNP. By default, nSPT=2, which means a hetSNP in a gamete can be imputed only if imputations are supported by more than 2 gametes and no imputation conflict from different supporting gametes. allowNA is the maximum number of gametes with NA values at a locus that are allowed after imputation. This parameter is useful only when the missing genotype rate is extremely high (eg., >80%) thus few markers can be completely imputed acroass all gamete cells. By default, allowNA=0. The markers, usually of a small number, with missing data that cannot be fully resolved by imputation are eliminated from the framework.

### imputation
imputedFrame <- hapiImupte(gmt = gmtFrame, nSPT = 2, allowNA = 0)

## Number of hetSNPs after imputation: 20994

head(imputedFrame)

##         gmt1 gmt2 gmt3 gmt4 gmt5
## 1135959    0    1    1    0    0
## 1136439    0    1    1    0    0
## 1136733    0    1    1    0    0
## 1137140    0    1    1    0    0
## 1137167    0    1    1    0    0
## 1137267    0    1    1    0    0

3.2.2.1 Majority voting for draft haplotype inference

With the imputed genotype data, draft haplotypes can be constructed based upon majority voting strategy implemented in the hapiPhase() function. The majority voting strategy is used to infer draft haplotypes by counting the number of links between two adjacent markers. Because recombination frequency is very low, the appearance of links that represents true haplotypes (hap-links) should be much more frequent than those generated by crossovers (cv-links). Draft haplotypes can be determined through voting for the major links in sliding windows of 2 hetSNPs along the whole chromosome.

### majority voting
draftHap <- hapiPhase(gmt = imputedFrame) ###
head(draftHap)

##         hap haplink cvlink tlink ratio
## 1135959   0       5      0     5     0
## 1136439   0       5      0     5     0
## 1136733   0       5      0     5     0
## 1137140   0       5      0     5     0
## 1137167   0       5      0     5     0
## 1137267   0       5      0     5     0

### check positions with cv-links
draftHap[draftHap$cvlink>=1,]

##           hap haplink cvlink tlink ratio
## 2665203     0       4      1     5  0.25
## 4685329     0       4      1     5  0.25
## 11730915    0       4      1     5  0.25
## 29414096    0       4      1     5  0.25
## 29465722    0       4      1     5  0.25
## 35178660    0       4      1     5  0.25
## 35322900    0       4      1     5  0.25
## 133530288   0       4      1     5  0.25
## 141518185   0       4      1     5  0.25
## 144823941   0       4      1     5  0.25

3.2.2.2 Proofreading of draft haplotypes

The first step of proofreading is to filter out small regions (eg. < 1M) with multiple cv-links. Erroneous/ambiguous genotypes in complex genomic regions will largely affect the accuracy of imputation and thus impact the determination of hap-links by majority voting. If multiple cv-links are identified in a very small region, hetSNPs in this region will be eliminated from the frame and new draft haplotypes will be inferred again.

### identification of clusters of cv-links
cvCluster <- hapiCVCluster(draftHap = draftHap, cvlink = 2)
cvCluster

## [1] left  right
## <0 rows> (or 0-length row.names)

### determine hetSNPs in small regions involving multiple cv-links
filter <- c()
for (i in 1:nrow(cvCluster)) {
    filter <- c(filter, which (rownames(draftHap) >= cvCluster$left[i] & 
        rownames(draftHap) <= cvCluster$right[i]))
}

length(filter)

## [1] 0

### filter out hetSNPs in complex regions and infer new draft haplotypes
if (length(filter) > 0) {
    imputedFrame <- imputedFrame[-filter, ]
    draftHap <- hapiPhase(imputedFrame)
} 

In the second step of proofreading, a Maximum Parsimony of Recombination (MPR) strategy is adopted to proofread draft haplotypes at genomic positions where there are potential incorrect cv-links. By MPR, we believe that the most likely haplotype should be the one generating less number of crossovers across all gamete cells. The cvlink parameter is used to determine the positions for dividing draft haplotypes into blocks for proofreading. With 5 or more gametes, positions that harboring more than 2 cv-links across all the gametes will be used by default. With 3 or 4 gametes, probably all the positions with cv-links should be proofread.

finalDraft <- hapiBlockMPR(draftHap = draftHap, gmtFrame = gmtFrame, cvlink = 1)

## Size of blocks: 132,119,600,3957,53,1340,75,13412,555,264,487

## 2 blocks are removed !

## Number of crossovers given haplotype 1/2: 1/4
## Number of crossovers given haplotype 1/2: 1/4
## Number of crossovers given haplotype 1/2: 1/4

## Number of crossovers given haplotype 1/2: 0/5
## Number of crossovers given haplotype 1/2: 0/5

## Number of crossovers given haplotype 1/2: 1/4
## Number of crossovers given haplotype 1/2: 1/4
## Number of crossovers given haplotype 1/2: 1/4

head(finalDraft)

## 1135959 1136439 1136733 1137140 1137167 1137267 
##       0       0       0       0       0       0