This is a test version of Biostars. For the public version, visit https://www.biostars.org.
Splicing QTL plot

I have performed splicing QTL (sQTL) analysis using leafcutter and FastQTL to get the output:

Phenotype_id    Number_of_variants_tested   MLE_Shape1  MLE_shape2  Dummy   variant_id  Distance    nominal_p-value slope   first_permutation_p-value   Second_permutation_p-value
1:3913:3996:clu_7241_NA 46  0.779368    10.5604 475.195 snp_1_2321  -1593   2.55028E-06 -0.704681   0.00789921  0.00345113

Now I am interested to plot SNPs along with intron coordinates to generate something like a Figure 2A

enter image description here

so I have subset the above file to get intron coordinates and SNP position and now the structure of my input file is like this:

chr Intron_start    Intron_end  SNP_position    Pvalue
1   3913    3996    2321    0.00345113
1   3913    4001    4313    0.0116419
1   7447    7564    7464    0.0160019
1   7450    7564    7465    0.0348276

I am able to generate a plot using Intron start and snp position using this command:

Trans <- read.delim("combine_benj_1000_final_5", header=TRUE, sep="\t")
theme_set(theme_bw())  # pre-set the bw theme.
# Scatterplot
g <- ggplot(Trans, aes(SNP_position, Intron_start))
g + geom_point(aes(col=Pvalue)) +
  labs(y="Gene Position",
       x="SNP Position")

Generated Plot

plot But I am not sure how can I use Intron coordinates (start and end) to plot against SNP_position, as for an intron or gene I think we need to specify its start and end location to correlate it with SNP position.

Secondly, how can I plot based on chromosomes?

Here is Input trial dataset

Any help will be highly appreciated.

sqtl snp ggplots r
• 2,310 views • link
updated by zx8754 zx8754 • written by waqaskhokhar999 waqaskhokhar999

As per my understanding, I have performed the trans - eQTL for small RNA based on genes associated SNPs , where each chromosome has shown. I am providing you the code , hope you may find it helpful. 

`rm(list=ls())
setwd("H:/Data/sRNA eQTL data/eQTL smallRNA related/RNA_data_results/17_6_17_srna/8_7_newetql/all_final/final_final/")


nini<-read.table("allclusterfilegwas.txt",head=F,sep="\t")

pplo1<-nini[which(nini[,1]<1e-15),]
huang<-cbind(-log10(pplo1[,1]),pplo1[,2:3],abs(pplo1[,3]-pplo1[,2]))



gg1<-huang[which(huang[,4]<=5000000),]


gg2<-huang[which(huang[,4]>5000000),]




hho<-rgb(red=0, green=0, blue = (0:15)/15,alpha=1)
pp2<-hho[c(8,10,12,14,16)]
hhn<-rgb(red=(0:15)/15, green =0, blue = 0,alpha=1)
pp1<-hhn[c(8,10,12,14,16)]


rr1<-NULL
rr1[intersect(which(gg1[,1]>=15),which(gg1[,1]<20))]<-pp1[5]
rr1[intersect(which(gg1[,1]>=20),which(gg1[,1]<25))]<-pp1[4]
rr1[intersect(which(gg1[,1]>=25),which(gg1[,1]<30))]<-pp1[3]
rr1[intersect(which(gg1[,1]>=30),which(gg1[,1]<35))]<-pp1[2]
rr1[intersect(which(gg1[,1]>=35),which(gg1[,1]<40))]<-pp1[1]





rr2<-NULL
rr2[intersect(which(gg2[,1]>=15),which(gg2[,1]<20))]<-pp2[5]
rr2[intersect(which(gg2[,1]>=20),which(gg2[,1]<25))]<-pp2[4]
rr2[intersect(which(gg2[,1]>=25),which(gg2[,1]<30))]<-pp2[3]
rr2[intersect(which(gg2[,1]>=30),which(gg2[,1]<35))]<-pp2[2]
rr2[intersect(which(gg2[,1]>=35),which(gg2[,1]<40))]<-pp2[1]

bbb<-data.frame(cbind(gg1,rr1))
bbb1<-data.frame(cbind(gg2,rr2))

colnames(bbb)<-1:5
colnames(bbb1)<-1:5
rownames(bbb)<-NULL
rownames(bbb1)<-NULL

ooo<-rbind(bbb1,bbb)





ttt<-read.table("eachchr_pos.txt",head=F,sep="\t")

plot(ooo[,2],ooo[,3],type="n",xlim=c(0,ttt[11,1]),ylim=c(0,ttt[11,1]),xlab="",ylab="",axe=F,xaxs="i",yaxs="i")
points(ooo[,2],ooo[,3],pch=20,cex=0.6,col=ooo[,5])
box(col="gray")
abline(v=ttt[2:10,1],col="gray")
abline(h=ttt[2:10,1],col="gray")

##plot(ooo[,2],ooo[,3],type="n",xlim=c(0,ttt[11,1]),ylim=c(0,ttt[11,1]),xlab="",ylab="",axe=F,xaxs="i",yaxs="i")
niu<-read.table("srna_gene.txt",head=F,sep="\t")

niu1<-cbind(niu[,1:2],ttt[niu[,4],1]+niu[,6],niu[,7])
##par(new=T)

points(niu1[,3],rep(2000000000,nrow(niu1)),pch=15,cex=0.3,col="orange")
text(niu1[,3],rep(600000000,nrow(niu1)),niu1[,2],cex=0.4,adj=1)`
• 0 views • link
written by abdul.suboor123 • abdul.suboor123

Read about cut, it will make your code 50% less.

• 0 views • link
written by zx8754 zx8754

can you please also provide me input dataset?

• 0 views • link
written by waqaskhokhar999 waqaskhokhar999

0 answers

No answers yet.

Log in to answer this question.

More posts like this

No similar posts found.