colnames(snp_position) <- X_chr1_cm[,1]
rownames(snp_position) <- NULL
myGD_ <- GD.candidate[,-1]
myGD_ <- X_chr1[,-1]
myGD_ <- replace(myGD_, myGD_ == 0, "AA")
myGD_ <- replace(myGD_, myGD_ == 1, "AB")
myGD_ <- replace(myGD_, myGD_ == 2, "BB")
myGD_with_posi <- rbind(snp_position,myGD_)
dim(myGD_)
dim(snp_position)
myGD_[1:5,1:5]
snp_position[1:5,1:5]
snp_position[1:2,1:5]
#get SNPs on the chromosome one
X=myGD[,-1]
index1=myGM[,2]==2
X_chr1 = X[,index1]
X_chr1_cm = myGM[index1,]
X_chr1_cm$Position = X_chr1_cm$Position / 1000000
# simulate Y with 10 SNPs in first five chromosomes
index1to5=myGM[,2]<6
X1to5 = X[,index1to5]
taxa=myGD[,1]
set.seed(99164)
GD_1to5=cbind(taxa,X1to5)
mySim=GAPIT.Phenotype.Simulation(GD=GD_1to5,GM=myGM[index1to5,],h2=.75,NQTN=20, effectunit =.95,QTNDist="normal")
# change the data format for R/qtl
snp_position <- t(X_chr1_cm[,2:3])
colnames(snp_position) <- X_chr1_cm[,1]
rownames(snp_position) <- NULL
myGD_ <- replace(myGD_, myGD_ == 0, "AA")
myGD_ <- replace(myGD_, myGD_ == 1, "AB")
myGD_ <- replace(myGD_, myGD_ == 2, "BB")
myGD_with_posi <- rbind(snp_position,myGD_)
dim(myGD_)
dim(snp_position)
dim(X_chr1)
myGD_ <- X_chr1
myGD_ <- replace(myGD_, myGD_ == 0, "AA")
myGD_ <- replace(myGD_, myGD_ == 1, "AB")
myGD_ <- replace(myGD_, myGD_ == 2, "BB")
myGD_with_posi <- rbind(snp_position,myGD_)
dim(myGD_)
myGD_[1:5,1:5]
snp_position[1:2,1:5]
myY <- mySim$Y
myY <- myY[,-1]
myY <- c(NA,NA,myY)
data_all <- cbind(myY, myGD_with_posi)
write.csv(data_all,file = "test.csv", row.names=FALSE,na = "")
# R/qtl read the data and calculate the LOD
cross_my <- read.cross("csv", file = "test.csv", genotypes = c("AA","AB","BB"))
out.mr <- mqmscan(cross_my,model = "additive")
# QTL significance by permutation
results <- mqmpermutation(cross_my, scanfunction=mqmscan,n.cluster=2, n.perm=100, batchsize=25)
#get SNPs on the chromosome one
X=myGD[,-1]
index1=myGM[,2]==3
X_chr1 = X[,index1]
X_chr1_cm = myGM[index1,]
X_chr1_cm$Position = X_chr1_cm$Position / 1000000
# simulate Y with 10 SNPs in first five chromosomes
index1to5=myGM[,2]<6
X1to5 = X[,index1to5]
taxa=myGD[,1]
set.seed(99164)
GD_1to5=cbind(taxa,X1to5)
mySim=GAPIT.Phenotype.Simulation(GD=GD_1to5,GM=myGM[index1to5,],h2=.75,NQTN=20, effectunit =.95,QTNDist="normal")
# change the data format for R/qtl
snp_position <- t(X_chr1_cm[,2:3])
colnames(snp_position) <- X_chr1_cm[,1]
rownames(snp_position) <- NULL
myGD_ <- X_chr1
myGD_ <- replace(myGD_, myGD_ == 0, "AA")
myGD_ <- replace(myGD_, myGD_ == 1, "AB")
myGD_ <- replace(myGD_, myGD_ == 2, "BB")
myGD_with_posi <- rbind(snp_position,myGD_)
myY <- mySim$Y
myY <- myY[,-1]
myY <- c(NA,NA,myY)
data_all <- cbind(myY, myGD_with_posi)
write.csv(data_all,file = "test.csv", row.names=FALSE,na = "")
# R/qtl read the data and calculate the LOD
cross_my <- read.cross("csv", file = "test.csv", genotypes = c("AA","AB","BB"))
out.mr <- mqmscan(cross_my,model = "additive")
#get SNPs on the chromosome one
X=myGD[,-1]
index1=myGM[,2]==10
X_chr1 = X[,index1]
X_chr1_cm = myGM[index1,]
X_chr1_cm$Position = X_chr1_cm$Position / 1000000
# simulate Y with 10 SNPs in first five chromosomes
index1to5=myGM[,2]<6
X1to5 = X[,index1to5]
taxa=myGD[,1]
set.seed(99164)
GD_1to5=cbind(taxa,X1to5)
mySim=GAPIT.Phenotype.Simulation(GD=GD_1to5,GM=myGM[index1to5,],h2=.75,NQTN=20, effectunit =.95,QTNDist="normal")
# change the data format for R/qtl
snp_position <- t(X_chr1_cm[,2:3])
colnames(snp_position) <- X_chr1_cm[,1]
rownames(snp_position) <- NULL
myGD_ <- X_chr1
myGD_ <- replace(myGD_, myGD_ == 0, "AA")
myGD_ <- replace(myGD_, myGD_ == 1, "AB")
myGD_ <- replace(myGD_, myGD_ == 2, "BB")
myGD_with_posi <- rbind(snp_position,myGD_)
myY <- mySim$Y
myY <- myY[,-1]
myY <- c(NA,NA,myY)
data_all <- cbind(myY, myGD_with_posi)
write.csv(data_all,file = "test.csv", row.names=FALSE,na = "")
# R/qtl read the data and calculate the LOD
cross_my <- read.cross("csv", file = "test.csv", genotypes = c("AA","AB","BB"))
out.mr <- mqmscan(cross_my,model = "additive")
# R/qtl read the data and calculate the LOD
cross_my <- read.cross("csv", file = "test.csv", genotypes = c("AA","AB","BB"))
out.mr <- mqmscan(cross_my,model = "additive")
# R/qtl read the data and calculate the LOD
cross_my <- read.cross("csv", file = "test.csv", genotypes = c("AA","AB","BB"))
out.mr <- mqmscan(cross_my,model = "additive")
# R/qtl read the data and calculate the LOD
cross_my <- read.cross("csv", file = "test.csv", genotypes = c("AA","AB","BB"))
out.mr <- mqmscan(cross_my,model = "additive")
# R/qtl read the data and calculate the LOD
cross_my <- read.cross("csv", file = "test.csv", genotypes = c("AA","AB","BB"))
out.mr <- mqmscan(cross_my,model = "additive")
# R/qtl read the data and calculate the LOD
cross_my <- read.cross("csv", file = "test.csv", genotypes = c("AA","AB","BB"))
out.mr <- mqmscan(cross_my,model = "additive")
# R/qtl read the data and calculate the LOD
cross_my <- read.cross("csv", file = "test.csv", genotypes = c("AA","AB","BB"))
out.mr <- mqmscan(cross_my,model = "additive")
# R/qtl read the data and calculate the LOD
cross_my <- read.cross("csv", file = "test.csv", genotypes = c("AA","AB","BB"))
out.mr <- mqmscan(cross_my,model = "additive")
# R/qtl read the data and calculate the LOD
cross_my <- read.cross("csv", file = "test.csv", genotypes = c("AA","AB","BB"))
out.mr <- mqmscan(cross_my,model = "additive")
# R/qtl read the data and calculate the LOD
cross_my <- read.cross("csv", file = "test.csv", genotypes = c("AA","AB","BB"))
out.mr <- mqmscan(cross_my,model = "additive")
# R/qtl read the data and calculate the LOD
cross_my <- read.cross("csv", file = "test.csv", genotypes = c("AA","AB","BB"))
out.mr <- mqmscan(cross_my,model = "additive")
# R/qtl read the data and calculate the LOD
cross_my <- read.cross("csv", file = "test.csv", genotypes = c("AA","AB","BB"))
out.mr <- mqmscan(cross_my,model = "additive")
# R/qtl read the data and calculate the LOD
cross_my <- read.cross("csv", file = "test.csv", genotypes = c("AA","AB","BB"))
out.mr <- mqmscan(cross_my,model = "additive")
# R/qtl read the data and calculate the LOD
cross_my <- read.cross("csv", file = "test.csv", genotypes = c("AA","AB","BB"))
out.mr <- mqmscan(cross_my,model = "additive")
plot.map(pull.map(cross_my),est.map(cross_my))
mamscan
mqmscan
# R/qtl read the data and calculate the LOD
cross_my <- read.cross("csv", file = "test.csv", genotypes = c("AA","AB","BB"))
out.mr <- mqmscan(cross_my,model = "additive")
write.csv(data_all,file = "test.csv", row.names=FALSE,na = "")
# R/qtl read the data and calculate the LOD
cross_my <- read.cross("csv", file = "test.csv", genotypes = c("AA","AB","BB"))
out.mr <- mqmscan(cross_my,model = "additive")
# R/qtl read the data and calculate the LOD
cross_my <- read.cross("csv", file = "listeria.csv", genotypes = c("AA","AB","BB"))
out.mr <- mqmscan(cross_my,model = "additive")
out.mr <- mqmscan(cross_my)
data(map10)
simcross <- sim.cross(map10, type="f2", n.ind=100, missing.prob=0.02)
geno.image(simcross)
augmentedcross <- mqmaugment(simcross, minprob=1.0)
out.mr <- mqmscan(augmentedcross)
#get SNPs on the chromosome one
X=myGD[,-1]
index1=myGM[,2]==1
X_chr1 = X[,index1]
num_sample <- dim(X_chr1)[1]
num_snp <- dim(X_chr1)[2]
# Get the recombination rate of the SNPs against the 1st SNP
SNP_1 = matrix(rep(X_chr1[,1],num_snp-1),ncol = num_snp-1)
SNP_other = X_chr1[,-1]
SNP_count = SNP_1 + SNP_other
dim(SNP_count)
SNP_count[1:5,1:5]
# R/qtl read the data and calculate the LOD
cross_my <- read.cross("csv", file = "test.csv", genotypes = c("AA","AB","BB"))
out.mr <- mqmscan(cross_my，window.size = 50)
out.mr <- mqmscan(cross_my, window.size = 50)
out.mr <- mqmscan(cross_my, window.size = 50, estimate.map = TRUE)
X_chr1[1:5,1:5]
dim(X_chr1)
#get SNPs on the chromosome one
X=myGD[,-1]
index1=myGM[,2]==1
X_chr1 = X[,index1]
num_sample <- dim(X_chr1)[1]
num_snp <- dim(X_chr1)[2]
X_chr1[1:5,1:5]
dim(X_chr1)
# Get the recombination rate of the SNPs against the 1st SNP
SNP_pre = X_chr1[,-1]
SNP_after = X_chr1[,-dim(X_chr1)[2]]
# Get the recombination rate of the SNPs against the 1st SNP
SNP_after = X_chr1[,-1]
SNP_pre = X_chr1[,-dim(X_chr1)[2]]
SNP_pre[1:5,1:5]
dim(SNP_pre)
dim(SNP_after)
SNP_count = SNP_pre + SNP_after * 3
SNP_count[1:5,1:5]
SNP_count[SNP_count[,1] == c(1,3,5,7)]
SNP_count[,1][SNP_count[,1] == c(1,3,5,7)]
SNP_count[,1] == c(1,3,5,7)
SNP_count[,1] %in% c(1,3,5,7)
SNP_count[,1][SNP_count[,1] %in% c(1,3,5,7)]
#get SNPs on the chromosome one
X=myGD[,-1]
index1=myGM[,2]==1
X_chr1 = X[,index1]
num_sample <- dim(X_chr1)[1]
num_snp <- dim(X_chr1)[2]
# Get the recombination rate of the SNPs against the 1st SNP
SNP_after = X_chr1[,-1]
SNP_pre = X_chr1[,-num_snp]
SNP_count = SNP_pre + SNP_after * 3
recomb = c()
for (i in 1:dim(SNP_count)[2]){
ct_1 = length(SNP_count[,i][SNP_count[,i] %in% c(1,3,5,7)])
ct_2 = length(SNP_count[,i][SNP_count[,i] %in% c(0,8)])
ct_3 = length(SNP_count[,i][SNP_count[,i] %in% c(2,6)])
ct = min(ct_2, ct_3) + ct_1
recomb = c(recomb,ct/num_sample)
}
length(recomb)
recomb[1:5]
recomb > 0.5
recomb[recomb > 0.5]
exp(1)
for (i in 1:dim(SNP_count)[2]){
ct_1 = length(SNP_count[,i][SNP_count[,i] %in% c(1,3,5,7)])
ct_2 = length(SNP_count[,i][SNP_count[,i] %in% c(0,8)])
ct_3 = length(SNP_count[,i][SNP_count[,i] %in% c(2,6)])
# recombination rate
r = (min(ct_2, ct_3) + ct_1)*100/num_sample
# Haldane's map function
d = 50 * log( (1/(1-2*r)), base = exp(1))
recomb = c(recomb,d)
}
log(4,ba)
log(4,base = exp(1))
recomb
for (i in 1:dim(SNP_count)[2]){
ct_1 = length(SNP_count[,i][SNP_count[,i] %in% c(1,3,5,7)])
ct_2 = length(SNP_count[,i][SNP_count[,i] %in% c(0,8)])
ct_3 = length(SNP_count[,i][SNP_count[,i] %in% c(2,6)])
# recombination rate
r = (min(ct_2, ct_3) + ct_1)/num_sample * 100
# Haldane's map function
d = 50 * log( (1/(1-2*r)), base = exp(1))
recomb = c(recomb,d)
}
recomb = c()
for (i in 1:dim(SNP_count)[2]){
ct_1 = length(SNP_count[,i][SNP_count[,i] %in% c(1,3,5,7)])
ct_2 = length(SNP_count[,i][SNP_count[,i] %in% c(0,8)])
ct_3 = length(SNP_count[,i][SNP_count[,i] %in% c(2,6)])
# recombination rate
r = (min(ct_2, ct_3) + ct_1)/num_sample * 100
# Haldane's map function
d = 50 * log( (1/(1-2*r)), base = exp(1))
recomb = c(recomb,d)
}
recomb
recomb = c()
rs = c()
for (i in 1:dim(SNP_count)[2]){
ct_1 = length(SNP_count[,i][SNP_count[,i] %in% c(1,3,5,7)])
ct_2 = length(SNP_count[,i][SNP_count[,i] %in% c(0,8)])
ct_3 = length(SNP_count[,i][SNP_count[,i] %in% c(2,6)])
# recombination rate
r = (min(ct_2, ct_3) + ct_1)/num_sample * 100
rs = c(rs,r)
# Haldane's map function
d = 50 * log( (1/(1-2*r)), base = exp(1))
recomb = c(recomb,d)
}
r
rs
50 * log( (1/(1-2*rs[1])), base = exp(1))
recomb = c()
rs = c()
for (i in 1:dim(SNP_count)[2]){
ct_1 = length(SNP_count[,i][SNP_count[,i] %in% c(1,3,5,7)])
ct_2 = length(SNP_count[,i][SNP_count[,i] %in% c(0,8)])
ct_3 = length(SNP_count[,i][SNP_count[,i] %in% c(2,6)])
# recombination rate
r = (min(ct_2, ct_3) + ct_1)/num_sample
rs = c(rs,r)
# Haldane's map function
d = 50 * log( (1/(1-2*r)), base = exp(1))
recomb = c(recomb,d)
}
recomb
50 * log( (1/(1-2*rs[1])), base = exp(1))
50 * log( (1/(1-2*rs[2])), base = exp(1))
1/(1-2*rs[2])
rs[2]
SNP_count[1:5,1:5]
(14+133)/281
50 * log( 0, base = exp(1))
50 * log( (1/(1-2*0.49)), base = exp(1))
50 * log( (1/(1-2*0.499)), base = exp(1))
50 * log( (1/(1-2*0.49999)), base = exp(1))
#Import data
myGD=read.table(file="http://zzlab.net/GAPIT/data/mdp_numeric.txt",head=T)
X.raw=myGD[,-1] #keep as original for comparison
X=X.raw # Create a new variable than may be changed later
myGD[1:5,1:5]
X[1:5,1:5]
#Set missing values
n=nrow(X)
m=ncol(X)
dp=m*n #total data points
uv=runif(dp)
hist(uv)
mr=.2 #missing rate
missing=uv<mr
length(missing)
length(uv)
missing[1:10]
#Format indicator as matrix
index.m=matrix(missing,n,m)
dim(index.m)
index.m[1:5,1:5]
#Set missing values as NA
X[index.m]=NA
X.raw[1:5,1:5]
X[1:5,1:5]
# X is the genotype data with n samples and m markers
n=nrow(X)
m=ncol(X)
fn=colSums(X, na.rm=T) # sum of genotypes for all individuals
dim(fn)
length(fn)
dim(X)
fn[1]
fc=colSums(floor(X/3+1),na.rm=T) #count number of non missing individuals
fc[1]
x[1:5,1]
X[1:5,1]
sum(X[1:5,1])
sum(X[1:5,1], na.rm = T)
sum(floor(X[1:5,1]/3+1), na.rm=T)
sum(X[1:5,1], na.rm = T) / (2*sum(floor(X[1:5,1]/3+1), na.rm=T))
floor(1/3)
floor(3/3)
floor(2/3)
floor(2/3 + 1)
hist(X[,1])
floor((c(1,1,NA)/3 + 1), na.rm=T)
floor((c(1,1,NA)/3 + 1))
floor(X/3+1)
a <- floor(X/3+1)
a[1:5,1:5]
fn[1]/(2*fc[1])
length(X[X[,1] == 2])
length(X[X[,1] == 2],1)
length(X[X[,1] == 2, 1])
220/length(X[,1])
hist(X[X[,1] == 2, 1])
X[X[,1] == 2, 1]
sum(c(0,1,2))
X[,1][1:10]
?runif
runif(nrow(X))
fn
fa
fc
fa=fn/(2*fc) #Frequency of allele "2"
fa[1]
fn[1]
fc[1]
338/(2*226)
index.a=runif(n)<fa[i]
index.a=runif(n)<fa[1]
index.a
fa[1]
is.na(X[,i])
is.na(X[,1])
X[1:5,1]
index.m2=index.a  &  index.na
index.na=is.na(X[,1])
index.m2=index.a  &  index.na
length(index.m2)
index.m2
index.m0=!index.a  &  index.na
index.m0
index.a
index.m2[1:5]
index.m0[1:5]
hist(index.m0)
index.m0
index.m2[1:10]
index.m0[1:10]
index.a[1:5]
index.a[1:3]
FALSE & TRUE
TRUE & TRUE
dim(index.m)
#Impute
XI= StochasticImpute(X)
return(X)}
#Define Stochastic Impute function
StochasticImpute = function(X){
# X is the genotype data with n samples and m markers
n=nrow(X)
m=ncol(X)
# Calculatet the frequency of allele A
fn=colSums(X, na.rm=T) # sum of genotypes for all individuals to get the number of A.
fc=colSums(floor(X/3+1),na.rm=T) #count number of non missing individuals
fa=fn/(2*fc) #Frequency of allele A.
for(i in 1:m){
index.na=is.na(X[,i])  # where is the missing value
index.a=runif(n)<fa[i] # randomly sign positions where is less than the frequency of allele "2"
index.m2 = index.na & index.a  # For those missing value, if it is less than the frequency of allele "2"
X[index.m2,i]=2                # then, sign the 2
index.m0 = index.na & !index.a # For those missing value, if it is not less than the frequency of allele "2"
X[index.m0,i]=0                # then, sign the 0
}
return(X)}
#Impute
XI= StochasticImpute(X)
index.m
index.m[1:5,1:5]
X.raw[index.m]
length(X.raw[index.m])
dp
173844/dp
nrep=100
myimp=replicate(nrep,{
uv=runif(dp)
#hist(uv)
missing=uv<mr
length(missing)
missing[1:10]
index.m=matrix(missing,n,m)
dim(index.m)
X[index.m]=NA
X.raw[1:5,1:5]
X[1:5,1:5]
#=======================================
#Impute with StochasticImpute
XI= StochasticImpute(X)
#Calcuate accuracy
accuracy.r=cor(X.raw[index.m], XI[index.m])
index.match=X.raw==XI
index.mm=index.match&index.m
accuracy.m=length(X[index.mm])/length(X[index.m])
accuracy.r
accuracy.m
acc=c(accuracy.r, accuracy.m)
})
dim(X)
#Impute using marker as features
X.knn= impute.knn(as.matrix((X)), k=10)
library(impute)
#Impute using marker as features
X.knn= impute.knn(as.matrix((X)), k=10)
plot(X.knn$data[index.m],pch=".")
str(X.knn)
X[1:5,1:5]
hist(X.knn$data[index.m])
index.m[1:5,1:5]
length(data[index.m])
length(X.knn$data[index.m])
X.knn$data[1:5,1:5]
#change to integer
X.knn.I= X.knn$data
Index1= X.knn$data<1
Index2= X.knn$data>1
X.knn.I[Index1]= 0
X.knn.I[Index2]= 2
#Correlation
accuracy.r=cor(X.raw[index.m], (X.knn.I)[index.m])
#Match proportion
index.match=X.raw== X.knn.I
index.mm=index.match&index.m
accuracy.m=length(X[index.mm])/length(X[index.m])
accuracy.r
accuracy.m
getwd()
#Convert to BEAGLE input format
index0=X==0
index1=X==1
index2=X==2
indexna=is.na(X)
X2=X
X2[index0]="A\tA"
X2[index1]="A\tB"
X2[index2]="B\tB"
X2[indexna]="?\t?"
myGD2=cbind("M",myGD[,1],X2)
myGD2[1:5,1:5]
X[1:3,1:3]
X[1:5,1:3]
setwd("/Users/ZhouTang/Documents/WSU_courses/2021 spring/545/2021/Lab/Lab4/lab4_imputation")
write.table(myGD2,file="test.bgl",quote=F,sep="\t",col.name=F,row.name=F)
#Impute with BEAGLE
system("java -Xmx12g -jar /Users/ZhouTang/Documents/WSU_courses/2021 spring/545/2021/Lab/Lab4/lab4_imputation/beagle.jar unphased=test.bgl missing=? out=test1" )
#Impute with BEAGLE
system("java -Xmx12g -jar /Users/ZhouTang/Documents/WSU_courses/2021\ spring/545/2021/Lab/Lab4/lab4_imputation/beagle.jar unphased=test.bgl missing=? out=test1" )
