setwd("/Users/matteo/Desktop/") library(dplyr) library(dunn.test) library(IRanges) library(GenomicRanges) ##### Pst Pst<-read.table("Pst_results_h2_0.5.txt",header=TRUE) #Read gene expression categories all2<-read.table("Allgenes_allclassifications.txt",header=TRUE) all2<-all2[,c(1,5)] all2$status <- as.factor(all2$status) genome4<-inner_join(Pst,all2,by="gene_id") #Kruskal-Wallis test kruskal.test(pst ~ status, data = genome4) #Dunn test dunn.test(x=genome4$pst,g=genome4$status,method="bonferroni",kw=TRUE) ##### Adxmiture proportions (fd) #read gene annotation (in chromosomal coordinates) annotation<-read.table("chromo_gene_info.txt",header=TRUE) withposition2<-inner_join(all2,annotation,by="gene_id") #Add admixture proportion estimates (100kb windows) admixture2<-read.csv("bar92.DP8MP4BIMAC2HET75.fourPopPol_melG_melW_cyd_num.w100m1s20.merged.csv",header=TRUE) admixture2<- subset(admixture2, start %% 100000 ==1) #non overlapping windows admixture2$fd[admixture2$fd < 0] <- 0 #negative values set to 0 #Merge admixture proportions and gene coordinates info GRgenes2<-makeGRangesFromDataFrame(withposition2, seqnames.field=c("scaffold")) #intermediate GRadmix<-makeGRangesFromDataFrame(admixture2, seqnames.field=c("scaffold")) olaps2 <- findOverlaps(GRgenes2, GRadmix) a2<-cbind(withposition2[queryHits(olaps2),], admixture2[subjectHits(olaps2),]) genome_INT<-a2 genome_INT<-genome_INT[,c(1:5,17)] #Dunn test dunn.test(x=genome_INT$fd,g=genome_INT$status,method="bonferroni",kw=TRUE)