plotHistBV <- function(nameTableGEBV, trait) {
# Plotting
dataGEBV = fread(paste0(nameTableGEBV, ".txt"))
sizeBin = (range(dataGEBV$GEBV) %>% diff())/20
plotGEBV = ggplot(data = dataGEBV, aes(x = GEBV, fill = ..density..)) +
labs(title = sprintf("Distribution of GEBVs for %s", trait)) +
geom_histogram(binwidth = sizeBin, colour = "black") +
theme(text = element_text(size = 20))
ggsave(plot = plotGEBV, file = paste0(nameTableGEBV, ".png"), height = 10, width = 14)
}
plotCV <- function(isRaw, nameTableRaw, nameTableAcc) {
# Plotting
if (isRaw) {
dataRaw = fread(paste0(nameTableRaw, ".txt"))
R2 = cor(dataRaw$obs, dataRaw$pre)^2
plotRaw = ggplot(data = dataRaw, aes(x = obs, y = pre)) +
stat_smooth(method = 'lm') +
geom_point() +
labs(title = sprintf("Observed Values Against Predictions Over Iterations (r2 = %.2f)", R2),
x = "Observed Values in Testing Dataset",
y = "Predicted Values in Testing Dataset") +
theme(text = element_text(size = 20))
ggsave(plot = plotRaw, file = paste0(nameTableRaw, ".png"), height = 10, width = 11)
}
dataAcc = fread(paste0(nameTableAcc, ".txt"))
dataAcc2 = data.table(value = c(dataAcc$r, dataAcc$rmse),
type = rep(c("Correlation (r)", "RMSE"), each = nrow(dataAcc)))
plotAcc = ggplot(data = dataAcc2, aes(x = value, fill = type)) +
labs(title = "Distribution of Prediction Accuracy") +
geom_density(alpha = .5) +
facet_grid(. ~ type, scales = "free") +
theme(text = element_text(size = 20))
ggsave(plot = plotAcc, file = paste0(nameTableAcc, ".png"), height = 10, width = 14)
}
`iPat.QQ` <-
function(P.values, plot.type = "log_P_values", filename = "_", name.of.trait = "Trait", DPP = 50000, plot.style = "rainbow"){
#Object: Make a QQ-Plot of the P-values
#Options for plot.type = "log_P_values" and "P_values"
#Output: A pdf of the QQ-plot
#Authors: Alex Lipka and Zhiwu Zhang
# Last update: May 9, 2011
##############################################################################################
# Sort the data by the raw P-values
P.values=P.values[!is.na(P.values)]
P.values=P.values[P.values>0]
P.values=P.values[P.values<=1]
if(length(P.values[P.values>0])<1) return(NULL)
N=length(P.values)
DPP=round(DPP/4) #Reduce to 1/4 for QQ plot
P.values <- P.values[order(P.values)]
#Set up the p-value quantiles
p_value_quantiles <- (1:length(P.values))/(length(P.values)+1)
if(plot.type == "log_P_values")
{
log.P.values <- -log10(P.values)
log.Quantiles <- -log10(p_value_quantiles)
index=GAPIT.Pruning(log.P.values,DPP=DPP)
log.P.values=log.P.values[index ]
log.Quantiles=log.Quantiles[index]
if(plot.style=="FarmCPU"){
pdf(paste0(filename, "_QQ-Plot.pdf"),width = 5,height=5)
par(mar = c(5,6,5,3))
}
if(plot.style=="rainbow"){
pdf(paste0(filename, "_QQ-Plot.pdf"), width = 5,height=5)
par(mar = c(5,6,5,3))
}
#Add conficence interval
N1=length(log.Quantiles)
## create the confidence intervals
c95 <- rep(NA,N1)
c05 <- rep(NA,N1)
for(j in 1:N1){
i=ceiling((10^-log.Quantiles[j])*N)
if(i==0)i=1
c95[j] <- qbeta(0.95,i,N-i+1)
c05[j] <- qbeta(0.05,i,N-i+1)
#print(c(j,i,c95[j],c05[j]))
}
#CI Lines
#plot(log.Quantiles, -log10(c05), xlim = c(0,max(log.Quantiles)), ylim = c(0,max(log.P.values)), type="l",lty=5, axes=FALSE, xlab="", ylab="",col="black")
#par(new=T)
#plot(log.Quantiles, -log10(c95), xlim = c(0,max(log.Quantiles)), ylim = c(0,max(log.P.values)), type="l",lty=5, axes=FALSE, xlab="", ylab="",col="black")
#CI shade
plot(NULL, xlim = c(0,max(log.Quantiles)), ylim = c(0,max(log.P.values)), type="l",lty=5, lwd = 2, axes=FALSE, xlab="", ylab="",col="gray")
index=length(c95):1
polygon(c(log.Quantiles[index],log.Quantiles),c(-log10(c05)[index],-log10(c95)),col='gray',border=NA)
#Diagonal line
abline(a = 0, b = 1, col = "red",lwd=2)
#data
par(new=T)
if(plot.style=="FarmCPU"){
plot(log.Quantiles, log.P.values, cex.axis=1.1, cex.lab=1.3, lty = 1, lwd = 2, col = "Black" ,bty='l', xlab =expression(Expected~~-log[10](italic(p))), ylab = expression(Observed~~-log[10](italic(p))), main = paste(name.of.trait,sep=""),pch=20)
}
if(plot.style=="rainbow"){
plot(log.Quantiles, log.P.values, xlim = c(0,max(log.Quantiles)), ylim = c(0,max(log.P.values)), cex.axis=1.1, cex.lab=1.3, lty = 1, lwd = 2, col = "Blue" ,xlab =expression(Expected~~-log[10](italic(p))),ylab = expression(Observed~~-log[10](italic(p))), main = paste(name.of.trait,sep=""))
}
dev.off()
}
if(plot.type == "P_values")
{
pdf(paste("QQ-Plot_", name.of.trait,".pdf" ,sep = ""))
par(mar = c(5,5,5,5))
qqplot(p_value_quantiles, P.values, xlim = c(0,1),
ylim = c(0,1), type = "l" , xlab = "Uniform[0,1] Theoretical Quantiles",
lty = 1, lwd = 1, ylab = "Quantiles of P-values from GWAS", col = "Blue",
main = paste(name.of.trait,sep=" "))
abline(a = 0, b = 1, col = "red")
dev.off()
}
#print("GAPIT.QQ accomplished successfully!")
}
`iPat.Manhattan` <-
function(GI.MP = NULL,GD=NULL, filename = "_", name.of.trait = "Trait",plot.type = "Genomewise",
DPP=50000,cutOff=0.01,band=5,seqQTN=NULL,plot.style="Oceanic",CG=NULL,plot.bin=10^9){
#Object: Make a Manhattan Plot
#Options for plot.type = "Separate_Graph_for_Each_Chromosome" and "Same_Graph_for_Each_Chromosome"
#Output: A pdf of the Manhattan Plot
#Authors: Alex Lipka, Zhiwu Zhang, Meng Li and Jiabo Wang
# Last update: Oct 10, 2016
#Add r2 between candidata SNP and other markers in on choromosome
##############################################################################################
#print("Manhattan ploting...")
#do nothing if null input
if(is.null(GI.MP)) return
#if(is.null(GD)) return
#print("Dimension of GI.MP")
#print(dim(GI.MP))
#print(head(GI.MP))
#print(tail(GI.MP))
#print(CG)
#seqQTN=c(300,1000,2500)
#Handler of lable paosition only indicated by negatie position
position.only=F
if(!is.null(seqQTN)){
if(seqQTN[1]<0){
seqQTN=-seqQTN
position.only=T
}
}
#if(is.null(GD)) print ("GD is not same dim as GM")
borrowSlot=4
GI.MP[,borrowSlot]=0 #Inicial as 0
GI.MP[,5]=1:(nrow(GI.MP))
if(!is.null(seqQTN))GI.MP[seqQTN,borrowSlot]=1
GI.MP=matrix(as.numeric(as.matrix(GI.MP) ) ,nrow(GI.MP),ncol(GI.MP))
GI.MP=GI.MP[order(GI.MP[,2]),]
GI.MP=GI.MP[order(GI.MP[,1]),]
#print("@@@@@")
#print(dim(GD))
#print(dim(GI.MP))
if(!is.null(GD))
{ if(ncol(GD)!=nrow(GI.MP))
{print("GD does not mach GM in Manhattan !!!")
return
}}
#print("!!")
#GI.MP[,5]=1:(nrow(GI.MP))
#print(head(GI.MP,20))
#Remove all SNPs that do not have a choromosome, bp position and p value(NA)
GI.MP <- GI.MP[!is.na(GI.MP[,1]),]
GI.MP <- GI.MP[!is.na(GI.MP[,2]),]
GI.MP <- GI.MP[!is.na(GI.MP[,3]),]
#Retain SNPs that have P values between 0 and 1 (not na etc)
GI.MP <- GI.MP[GI.MP[,3]>0,]
GI.MP <- GI.MP[GI.MP[,3]<=1,]
if(!is.null(GD)) GD=GD[,GI.MP[,3]<=1]
#Remove chr 0 and 99
GI.MP <- GI.MP[GI.MP[,1]!=0,]
#GI.MP <- GI.MP[GI.MP[,1]!=99,]
#print(dim(GI.MP))
#print("Dimension of GI.MP after QC")
#print(dim(GI.MP))
numMarker=nrow(GI.MP)
bonferroniCutOff=-log10(cutOff/numMarker)
#Replace P the -log10 of the P-values
if(!is.null(GD))
{ if(ncol(GD)!=nrow(GI.MP))
{print("GD does not mach GM in Manhattan !!!")
return
}}
GI.MP[,3] <- -log10(GI.MP[,3])
index_GI=GI.MP[,3]>0
GI.MP <- GI.MP[index_GI,]
if(!is.null(GD)) GD=GD[,index_GI]
GI.MP[,5]=1:(nrow(GI.MP))
y.lim <- ceiling(max(GI.MP[,3]))
chm.to.analyze <- unique(GI.MP[,1])
#print(dim(GI.MP))
#print(dim(GD))
#print("name of chromosomes:")
#print(chm.to.analyze)
chm.to.analyze=chm.to.analyze[order(chm.to.analyze)]
numCHR= length(chm.to.analyze)
#GI.MP[,5]=1:(nrow(GI.MP))
bin.mp=GI.MP[,1:3]
bin.mp[,3]=0 # for r2
bin.mp[,1]=as.numeric(as.vector(GI.MP[,2]))+as.numeric(as.vector(GI.MP[,1]))*(10^(max(GI.MP[,1])+1))
#as.numeric(as.vector(GP[,3]))+as.numeric(as.vector(GP[,2]))*MaxBP
bin.mp[,2]=floor(bin.mp[,1]/plot.bin)
if(!is.null(GD)) X=GD
#print(head(bin.mp))
#Chromosomewise plot
if(plot.type == "Chromosomewise"&!is.null(GD))
{
#print("Manhattan ploting Chromosomewise")
GI.MP=cbind(GI.MP,bin.mp)
pdf(paste0(filename,"_Manhattan.Plot.Chromosomewise.pdf"), width = 10)
#par(mar = c(5,5,4,3), lab = c(8,5,7))
layout(matrix(c(1,1,2,1,1,1,1,1,1),3,3,byrow=TRUE), c(2,1), c(1,1), TRUE)
for(i in 1:numCHR)
{
#Extract SBP on this chromosome
subset=GI.MP[GI.MP[,1]==chm.to.analyze[i],]
# print(head(subset))
subset[,1]=1:(nrow(subset))
#sub.bin.mp=bin.mp[GI.MP[,1]==chm.to.analyze[i],]
#subset=cbind(subset,sub.bin.mp)
sig.mp=subset[subset[,3]>bonferroniCutOff,]
sig.index=subset[,3]>bonferroniCutOff ### index of significont SNP
num.row=nrow(sig.mp)
if(!is.null(dim(sig.mp)))sig.mp=sig.mp[!duplicated(sig.mp[,7]),]
num.row=nrow(sig.mp)
if(is.null(dim(sig.mp))) num.row=1
bin.set=NULL
r2_color=matrix(0,nrow(subset),2)
#r2_color
print(paste("select ",num.row," candidate gene in ",i," chromosome ",sep="") )
#print(sig.mp)
if(length(unique(sig.index))==2)
{
for(j in 1:num.row)
{ sig.mp=matrix(sig.mp,num.row,8)
#print(sig.mp[j,7])
#print(unique(subset[,7]))
bin.store=subset[which(subset[,7]==sig.mp[j,7]),]
if(is.null(dim(bin.store)))
{subset[which(subset[,7]==sig.mp[j,7]),8]=1
next
}
bin.index=unique(bin.store[,5])
subGD=X[,bin.store[,5]]
#print(dim(bin.store))
if(is.null(CG))candidata=bin.store[bin.store[,3]==max(bin.store[,3]),5]
if(length(candidata)!=1)candidata=candidata[1]
for (k in 1:ncol(subGD))
{
r2=cor(X[,candidata],subGD[,k])^2
#print(r2)
bin.store[k,8]=r2
}
#print(bin.store)
#r2_storage[is.na(r2_storage)]=0
#print(bin.store)
subset[bin.store[,1],8]=bin.store[,8]
#print()
}###end for each sig.mp
#sub.bin.mp=bin.mp[subset[,3]>bonferroniCutOff,]
#print(head(bin.set))
}###end if empty of sig.mp
#print("@@@@@@@@@@@@@@@@")
rm(sig.mp,num.row)
#print(head(subset))
#print(head(subset))
#print(dim(X))
y.lim <- ceiling(max(subset[,3]))+1 #set upper for each chr
if(length(subset)>3){
x <- as.numeric(subset[,2])/10^(6)
y <- as.numeric(subset[,3])
}else{
x <- as.numeric(subset[2])/10^(6)
y <- as.numeric(subset[3])
}
##print(paste("befor prune: chr: ",i, "length: ",length(x),"max p",max(y), "min p",min(y), "max x",max(x), "Min x",min(x)))
n_col=10
r2_color[,2]=subset[,8]
do_color=colorRampPalette(c("orangeRed", "blue"))(n_col)
#Prune most non important SNPs off the plots
order=order(y,decreasing = TRUE)
y=y[order]
x=x[order]
r2_color=r2_color[order,]
index=GAPIT.Pruning(y,DPP=round(DPP/numCHR))
x=x[index]
y=y[index]
r2_color=r2_color[index,]
r2_color[which(r2_color[,2]<=0.2),2]=do_color[n_col]
r2_color[which(r2_color[,2]<=0.4&r2_color[,2]>0.2),2]=do_color[n_col*0.8]
r2_color[which(r2_color[,2]<=0.6&r2_color[,2]>0.4),2]=do_color[n_col*0.6]
r2_color[which(r2_color[,2]<=0.8&r2_color[,2]>0.6),2]=do_color[n_col*0.4]
r2_color[which(r2_color[,2]<=1&r2_color[,2]>0.8),2]=do_color[n_col/n_col]
#print(unique(r2_color[,2]))
##print(paste("after prune: chr: ",i, "length: ",length(x),"max p",max(y), "min p",min(y), "max x",max(x), "Min x",min(x)))
par(mar=c(0,0,0,0))
par(mar=c(5,5,2,1),cex=0.8)
plot(y~x,type="p", ylim=c(0,y.lim), xlim = c(min(x), max(x)),
col = r2_color[,2], xlab = expression(Base~Pairs~(x10^-6)),
ylab = "-Log Base 10 p-value", main = paste("Chromosome",chm.to.analyze[i],sep=" "),
cex.lab=1.6,pch=21,bg=r2_color[,2])
abline(h=bonferroniCutOff,col="forestgreen")
##print("manhattan plot (chr) finished")
#layout.show(nf)
#provcol <-c("darkblue","cyan","green3","brown1","brown1")
#provcol <-heat.colors(50)
#par(mar=c(0,0,0,0))
par(mar=c(15,5,6,5),cex=0.5)
barplot(matrix(rep(1,times=n_col),n_col,1),beside=T,col=do_color,border=do_color,axes=FALSE,)
#legend(x=10,y=2,legend=expression(R^"2"),,lty=0,cex=1.3,bty="n",bg=par("bg"))
axis(3,seq(11,1,by=-2),seq(0,1,by=0.2))
}# end plot.type == "Chromosomewise"&!is.null(GD)
dev.off()
print("manhattan plot on chromosome finished")
} #Chromosomewise plot
#Genomewise plot
if(plot.type == "Genomewise")
{
#print("Manhattan ploting Genomewise")
#Set corlos for chromosomes
#nchr=max(chm.to.analyze)
nchr=length(chm.to.analyze)
#Set color schem
ncycle=ceiling(nchr/band)
ncolor=band*ncycle
#palette(rainbow(ncolor+1))
cycle1=seq(1,nchr,by= ncycle)
thecolor=cycle1
for(i in 2:ncycle){thecolor=c(thecolor,cycle1+(i-1))}
col.Rainbow=rainbow(ncolor+1)[thecolor]
col.FarmCPU=rep(c("#CC6600","deepskyblue","orange","forestgreen","indianred3"),ceiling(numCHR/5))
col.Rushville=rep(c("orangered","navyblue"),ceiling(numCHR/2))
col.Congress=rep(c("deepskyblue3","firebrick"),ceiling(numCHR/2))
col.Ocean=rep(c("steelblue4","cyan3"),ceiling(numCHR/2))
col.PLINK=rep(c("gray10","gray70"),ceiling(numCHR/2))
col.Beach=rep(c("turquoise4","indianred3","darkolivegreen3","red","aquamarine3","darkgoldenrod"),ceiling(numCHR/5))
#col.Oceanic=rep(c( '#EC5f67', '#F99157', '#FAC863', '#99C794', '#5FB3B3', '#6699CC', '#C594C5', '#AB7967'),ceiling(numCHR/8))
#col.Oceanic=rep(c( '#EC5f67', '#FAC863', '#99C794', '#6699CC', '#C594C5', '#AB7967'),ceiling(numCHR/6))
col.Oceanic=rep(c( '#EC5f67', '#FAC863', '#99C794', '#6699CC', '#C594C5'),ceiling(numCHR/5))
col.cougars=rep(c( '#990000', 'dimgray'),ceiling(numCHR/2))
if(plot.style=="Rainbow")plot.color= col.Rainbow
if(plot.style =="FarmCPU")plot.color= col.Rainbow
if(plot.style =="Rushville")plot.color= col.Rushville
if(plot.style =="Congress")plot.color= col.Congress
if(plot.style =="Ocean")plot.color= col.Ocean
if(plot.style =="PLINK")plot.color= col.PLINK
if(plot.style =="Beach")plot.color= col.Beach
if(plot.style =="Oceanic")plot.color= col.Oceanic
if(plot.style =="cougars")plot.color= col.cougars
#FarmCPU uses filled dots
mypch=1
if(plot.style =="FarmCPU")mypch=20
GI.MP <- GI.MP[order(GI.MP[,2]),]
GI.MP <- GI.MP[order(GI.MP[,1]),]
ticks=NULL
lastbase=0
#print("Manhattan data sorted")
#print(chm.to.analyze)
#change base position to accumulatives (ticks)
for (i in chm.to.analyze)
{
index=(GI.MP[,1]==i)
ticks <- c(ticks, lastbase+mean(GI.MP[index,2]))
GI.MP[index,2]=GI.MP[index,2]+lastbase
lastbase=max(GI.MP[index,2])
}
#print("Manhattan chr processed")
#print(length(index))
#print(length(ticks))
#print((ticks))
#print((lastbase))
x0 <- as.numeric(GI.MP[,2])
y0 <- as.numeric(GI.MP[,3])
z0 <- as.numeric(GI.MP[,1])
position=order(y0,decreasing = TRUE)
index0=GAPIT.Pruning(y0[position],DPP=DPP)
index=position[index0]
x=x0[index]
y=y0[index]
z=z0[index]
#Extract QTN
QTN=GI.MP[which(GI.MP[,borrowSlot]==1),]
#Draw circles with same size and different thikness
size=1 #1
ratio=10 #5
base=1 #1
themax=ceiling(max(y))
themin=floor(min(y))
wd=((y-themin+base)/(themax-themin+base))*size*ratio
s=size-wd/ratio/2
#print("Manhattan XY created")
####xiaolei update on 2016/01/09
if(plot.style =="FarmCPU"){
pdf(paste0(filename,"_Manhattan.Plot.Genomewise.pdf"), width = 13, height=5.75)
}else{
pdf(paste0(filename,"_Manhattan.Plot.Genomewise.pdf"), width = 13, height=5.75)
}
par(mar = c(3,6,5,1))
plot(y~x,xlab="",ylab=expression(-log[10](italic(p))) ,
cex.axis=1.5, cex.lab=2, ,col=plot.color[z],axes=FALSE,type = "p",pch=mypch,lwd=wd,cex=s+.3,main = paste(name.of.trait,sep=" "),cex.main=2.5)
#Label QTN positions
if(is.vector(QTN)){
if(position.only){abline(v=QTN[2], lty = 2, lwd=1.5, col = "grey")}else{
points(QTN[2], QTN[3], type="p",pch=21, cex=2,lwd=1.5,col="dimgrey")
points(QTN[2], QTN[3], type="p",pch=20, cex=1,lwd=1.5,col="dimgrey")
}
}else{
if(position.only){abline(v=QTN[,2], lty = 2, lwd=1.5, col = "grey")}else{
points(QTN[,2], QTN[,3], type="p",pch=21, cex=2,lwd=1.5,col="dimgrey")
points(QTN[,2], QTN[,3], type="p",pch=20, cex=1,lwd=1.5,col="dimgrey")
}
}
#Add a horizontal line for bonferroniCutOff
abline(h=bonferroniCutOff,col="forestgreen")
#Set axises
axis(1, at=ticks,cex.axis=1.5,labels=chm.to.analyze,tick=F)
axis(2, at=1:themax,cex.axis=1.5,labels=1:themax,tick=F)
box()
palette("default")
dev.off()
#print("Manhattan done Genomewise")
} #Genomewise plot
print("GAPIT.Manhattan accomplished successfully!zw")
} #end of GAPIT.Manhattan
`iPat.Genotype.View` <-function(myGD = NULL, filename = "_"){
# Object: Analysis for Genotype data:Distribution of SNP density,Accumulation,Moving Average of density,result:a pdf of the scree plot
# myG:Genotype data
# chr: chromosome value
# w1_start:Moving Average windows Start Position
# w1_end:Moving Average windows End Position
# mav1:Moving Average set value length
# Authors: You Tang and Zhiwu Zhang
# Last update: March 11, 2016
##############################################################################################
#heterozygosity of individuals and SNPs (By Zhiwu Zhang)
#print("Heterozygosity of individuals and SNPs (By Zhiwu Zhang)")
X=myGD[,-1]
H=1-abs(X-1)
het.ind=apply(H,1,mean)
het.snp=apply(H,2,mean)
ylab.ind=paste("Frequency (out of ",length(het.ind)," individuals)",sep="")
ylab.snp=paste("Frequency (out of ",length(het.snp)," markers)",sep="")
pdf(paste0(filename, "_heterozygosity.pdf"), width = 10, height = 6)
par(mfrow=c(1,2),mar=c(5,5,1,1)+0.1)
hist(het.ind,col="gray", main="",ylab=ylab.ind, xlab="Heterozygosity of individuals")
hist(het.snp,col="gray", main="",ylab=ylab.snp, xlab="Heterozygosity of markers")
dev.off()
}
`iPat.Phenotype.View` <-function(myY = NULL, filename = "_"){
# Object: Analysis for Phenotype data:Distribution of density,Accumulation,result:a pdf of the scree plot
# myY:Phenotype data
# Authors: You Tang
# Last update: Sep 7, 2015
##############################################################################################
print("GAPIT.Phenotype.View in press...")
if(is.null(myY)){stop("Validation Invalid. Please select read valid Phenotype flies !")}
y<-myY[!is.na(myY[,2]),2]
obs<-as.matrix(y)
traitname=colnames(myY)[2]
pdf(paste0(filename, "_phnotype_view.pdf"), width = 10, height = 6)
par(mar = c(5,5,5,5))
par(mfrow=c(2,2))
plot(obs,pch=1)
#hist(obs)
hist(obs,xlab="Density",main="",breaks=12, cex.axis=1,col = "gray")
boxplot(obs)
plot(ecdf(obs),col="red",bg="lightgray",xlab="Density",ylab="Accumulation",main="")
dev.off()
}
`Blink.LDRemoveDivided`<-function(GDneo=NULL,LD=NULL,Porder=NULL,bound=FALSE,model="A",orientation=NULL,l=NULL){
#Objects: LD remove, especially length(Porder)>10000
#Authors: Yao Zhou
#Last update: 08/15/2016
seqQTN=NULL
lp=length(Porder)
k=ceiling(lp/l)
GDneo=as.matrix(GDneo)
if(min(ncol(GDneo),nrow(GDneo))<201) bound=FALSE
if(orientation=="col"){
n=nrow(GDneo)
if(bound){
GDneo=GDneo[sample(n,200,replace=F),]
}
}else{
n=ncol(GDneo)
if(bound){
GDneo=GDneo[,sample(n,200,replace=F)]
}
GDneo=t(GDneo)
}
for(i in 1:k){
bottom=(i-1)*l+1
up=l*i
if(up>lp) up = lp
Porderb=Porder[bottom:up]
GDneob = GDneo[,bottom:up]
seqQTNs=Blink.LDRemoveBlock(GDneo=GDneob,LD=LD,Porder=Porderb,orientation="col",model=model)
seqQTN=append(seqQTN,seqQTNs)
}
rm(GDneob,seqQTNs,Porderb)
return(seqQTN)
}
`Blink.LDRemoveBlock`<-function(GDneo=NULL,LD=NULL,Porder=NULL,bound=FALSE,model="A",orientation=NULL){
#`Blink.LDRemove`<-function(GDneo=NULL,LD=NULL,Porder=NULL,bound=FALSE,model="A",orientation=NULL){
#Objects: Calculate LD and remove the correlated SNPs
#Authors: Yao Zhou
#Last Update: 03/03/16
if (model=="D"){
GDneo=1-abs(GDneo-1)
}
GDneo=as.matrix(GDneo)
if(min(ncol(GDneo),nrow(GDneo))<201) bound=FALSE
if(orientation=="col"){
n=nrow(GDneo)
if(bound){
GDneo=GDneo[sample(n,200,replace=F),]
}
}else{
n=ncol(GDneo)
if(bound){
GDneo=GDneo[,sample(n,200,replace=F)]
}
GDneo=t(GDneo)
}
corr=cor(GDneo)
corr[is.na(corr)]=1
corr[abs(corr)<=LD]=0
corr[abs(corr)>LD]=1
Psort=as.numeric(matrix(1,1,ncol(corr)))
for(i in 2:ncol(corr)){
p.a=Psort[1:(i-1)]
p.b=as.numeric(corr[1:(i-1),i])
index=(p.a==p.b)
index[(p.a==0)&(p.b==0)]=FALSE
if(sum(index)!=0) Psort[i]=0
}
seqQTN=Porder[Psort==1]
return(seqQTN)
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.