Was this page helpful?

G-statistics and plot

    Table of contents
    No headers
     
    test <- read.csv ("/Users/jieyin/Desktop/F2BSA/2858SNP4gscore.csv", sep=",",header=T)
    head(test)
    dim(test)
    c1 <- test$n1
    c2 <- test$n2
    c3 <- test$n3
    c4 <- test$n4
    chr <-test$chr
    pos <- test$pos
    ref<- test$ref
    alt <- test$alt
    #calculate mean coverage
    c1mean<-mean(c1)
    c2mean<-mean(c2)
    c3mean<-mean(c3)
    c4mean<-mean(c4)
    c<-2*mean(c1mean,c2mean,c3mean,c4mean)
    c
    dim(test)
    gvec<-rep(0, dim(test)[1])
    pval<-rep(0, dim(test)[1])
    logvec<-rep(NA, dim(test)[1])
     
    for ( i in 1:dim(test)[1]) {
    n1 <- test[i,6]
        n2 <- test[i,7]
        n3 <- test[i,8]
        n4 <- test[i,9]
    n1exp <-(n1+n2)*(n1+n3)/(n1+n2+n3+n4)
        n2exp <-(n1+n2)*(n2+n4)/(n1+n2+n3+n4)
        n3exp <-(n3+n4)*(n1+n3)/(n1+n2+n3+n4)
        n4exp <-(n3+n4)*(n2+n4)/(n1+n2+n3+n4)
     
        g1<-log(n1/n1exp)
        g2<-log(n2/n2exp)
        g3<-log(n3/n3exp)
        g4<-log(n4/n4exp)
        g<- 2*(n1*g1+n2*g2+n3*g3+n4*g4)
        
        p<- (1- pchisq(g,1))
        logp<- -log10(p)
        gvec[i]<-g
        pval[i]<-p
       
    }
    #gvec
    #pval
    #use FDR to find out the threshold 
    pvals<-pval
    x<-pvals
    x[is.nan(x)]<-0
    pvals<-x
    padjust<-p.adjust(pval, "fdr", n=length(pval))
     
     
    ############################## plot g-score against chromosome postion ####################################
     
     
    Gstat <- read.table("/home/jyin/mydir/AllSNP/WUE_allSNP_10kb.csv",sep=",",header=T)
    pval<-Gstat$pval
     
     
    x <- Gstat$gscore
    y <- Gstat$chr
    z <- Gstat$pos
    u <- Gstat$idx
    v <-Gstat$posidx
     
     xv<-c(0,28704,51603,76503,101526,127294,156675,185975)
     xl<-c(0,28704,22899,24900,25023,25768,29381,29300)
     yv<-c(0,20,40,60,80,100,120)
     yl<-c(seq(from = 0, to = 120, by = 20))
     xv2<-c(14352,40153,64053,89015,114410,141985,171325)
     xl2<-c("chr1","chr2","chr3","chr4","chr5","chr6","chr1")
     
     axis(side = 1, at = xv, tck = 0.01,labels=xl)
     axis(side = 2, at = yv, tck=0.01,labels= yl) 
     
     
     mtext(xl2, side = 1, line = -1.0, at = xv2, font = 2)
     
     # Add y-axis (tick marks and labels)
     #mtext("", side = 2, line = -1.20, at = seq(from = 0, to = 120, by = 20), font = 2)
     
     
     
     plot(x[y == 1]~ v[y==1], type="l", col="black",frame=TRUE, axes=FALSE,ylab="G-statistics",xlab="Physical Map (Kb)",xlim=c(0,max(v)), ylim=c(0,120)) 
     
     points(x[y==2]~v[y==2], type="l", col="black")
     points(x[y==3]~v[y==3], type="l", col="black")
     points(x[y==4]~v[y==4], type="l", col="black")
     points(x[y==5]~v[y==5], type="l", col="black")
     points(x[y==6]~v[y==6], type="l", col="black")
     points(x[y==7]~v[y==7], type="l", col="black")
     
     xv<-c(0,28704,51603,76503,101526,127294,156675,185975)
     xl<-c(0,28704,22899,24900,25023,25768,29381,29300)
     yv<-c(0,20,40,60,80,100,120)
     yl<-c(seq(from = 0, to = 120, by = 20))
     xv2<-c(14352,40153,64053,89015,114410,141985,171325)
     xl2<-c("chr1","chr2","chr3","chr4","chr5","chr6","chr1")
     
     axis(side = 1, at = xv, tck = 0.01,labels=xl)
     axis(side = 2, at = yv, tck= 0.01,labels= yl) 
     
     
     mtext(xl2, side = 1, line = -1.0, at = xv2, font = 2)
     
     # Add y-axis (tick marks and labels)
     #mtext("", side = 2, line = -1.20, at = seq(from = 0, to = 120, by = 20), font = 2)
     
     
     
    abline(h=gthr,col='red')           
    axis(2,at=gthr, line=-2, hadj=-0.5, label=gthr)        
    Was this page helpful?
    Tag page (Edit tags)
    • No tags
    You must login to post a comment.