coverage

#read in the files

library(data.table) #much faster than normal read in functions

males<-fread("scaf1631_males_depth.txt")

females<-fread("scaf1631_females_depth.txt")

dim(males)

dim(females)

males<-males[,-1]

females<-females[,-1]

#take means across columns

mmean<-rowMeans(males[,-1])

fmean<-rowMeans(females[,-1])

mp_males<-cbind(males[,1],mmean)

mp_fem<-cbind(females[,1],fmean)

#take means across windows of positions

bin <- seq(1,21870361,10000)

males_winmeans<-sapply(1:(length(bin)-1),

function(i)

sum(mp_males[mp_males[,1]>=bin[i] &

mp_males[,1] < bin[i+1], "mmean"]))

#fbin<-seq(1,21870361,10000)

fem_winmeans<-sapply(1:(length(bin)-1),

function(i)

sum(mp_fem[mp_fem[,1]>=bin[i] &

mp_fem[,1] < bin[i+1], "fmean"]))

plot(males_winmeans, fem_winmeans, col=c("red","blue"))

plot(males_winmeans/fem_winmeans,type='l')

abline(h=1)

abline(h=2,col="blue")

################## scaffold 1628 ####################################

males28<-fread("scaf1628_males_depth.txt")

females28<-fread("scaf1628_females_depth.txt")

dim(males28) #5457508 481

dim(females28) #4558600 354

males28<-males28[,-1]

females28<-females28[,-1]

#take means across columns

mmean28<-rowMeans(males28[,-1])

fmean28<-rowMeans(females28[,-1])

mp_males28<-cbind(males28[,1],mmean28)

mp_fem28<-cbind(females28[,1],fmean28)

mp_males28<-as.data.frame(mp_males28)

mp_fem28<-as.data.frame(mp_fem28)

#get the max position distance

max(mp_males45$V2) #23101924

max(mp_fem45$V2) #23101837

#take means across windows of positions

bin28 <- seq(1,23101924,10000)

males_winmeans28<-sapply(1:(length(bin28)-1),

function(i)

sum(mp_males28[mp_males28[,1]>=bin28[i] &

mp_males28[,1] < bin28[i+1], "mmean28"]))

#fbin<-seq(1,21870361,10000)

fem_winmeans28<-sapply(1:(length(bin28)-1),

function(i)

sum(mp_fem28[mp_fem28[,1]>=bin28[i] &

mp_fem28[,1] < bin28[i+1], "fmean28"]))

plot(males_winmeans28, fem_winmeans28, col=c("red","blue"))

plot(males_winmeans28/fem_winmeans28,type='l')

abline(h=1)

abline(h=2,col="blue")

################# scaffold 1645 ####################################

males45<-fread("scaf1645_males_depth.txt")

females45<-fread("scaf1645_females_depth.txt")

dim(males45) #4064187 481

dim(females45) #3406437 354

males45<-males45[,-1]

females45<-females45[,-1]

#take means across columns

mmean45<-rowMeans(males45[,-1])

fmean45<-rowMeans(females45[,-1])

mp_males45<-cbind(males45[,1],mmean45)

mp_fem45<-cbind(females45[,1],fmean45)

mp_males45<-as.data.frame(mp_males45)

mp_fem45<-as.data.frame(mp_fem45)

#get the max position distance

max(mp_males45$V2) #18387850

max(mp_fem45$V2) #18387815

#take means across windows of positions

bin45 <- seq(1,18387850,10000)

males_winmeans45<-sapply(1:(length(bin45)-1),

function(i)

sum(mp_males45[mp_males45[,1]>=bin45[i] &

mp_males45[,1] < bin45[i+1], "mmean45"]))

#fbin<-seq(1,21870361,10000)

fem_winmeans45<-sapply(1:(length(bin45)-1),

function(i)

sum(mp_fem45[mp_fem45[,1]>=bin45[i] &

mp_fem45[,1] < bin45[i+1], "fmean45"]))

plot(males_winmeans45, fem_winmeans45, col=c("red","blue"))

plot(males_winmeans45/fem_winmeans45,type='l')

abline(h=1)

abline(h=2,col="blue")

################# scaffold 1646 ####################################

males46<-fread("scaf1646_males_depth.txt")

females46<-fread("scaf1646_females_depth.txt")

dim(males46) #5467508 481

dim(females46) #4658600 354

males46<-males46[,-1]

females46<-females46[,-1]

#take means across columns

mmean46<-rowMeans(males46[,-1])

fmean46<-rowMeans(females46[,-1])

mp_males46<-cbind(males46[,1],mmean46)

mp_fem46<-cbind(females46[,1],fmean46)

mp_males46<-as.data.frame(mp_males46)

mp_fem46<-as.data.frame(mp_fem46)

#get the max position distance

max(mp_males46$V2) #18436971

max(mp_fem46$V2) #18436929

#take means across windows of positions

bin46 <- seq(1,18436971,10000)

males_winmeans46<-sapply(1:(length(bin46)-1),

function(i)

sum(mp_males46[mp_males46[,1]>=bin46[i] &

mp_males46[,1] < bin46[i+1], "mmean46"]))

#fbin<-seq(1,21870361,10000)

fem_winmeans46<-sapply(1:(length(bin46)-1),

function(i)

sum(mp_fem46[mp_fem46[,1]>=bin46[i] &

mp_fem46[,1] < bin46[i+1], "fmean46"]))

plot(males_winmeans46, fem_winmeans46, col=c("red","blue"))

plot(males_winmeans46/fem_winmeans46,type='l')

abline(h=1)

abline(h=2,col="blue")

#### final plot #############

pdf("coverage_scaffolds.pdf",width=15,height=7,bg="white")

par(mfrow=c(2,2))

par(mar=c(5,7,3,1))

#plot1

plot(males_winmeans/fem_winmeans,type='l', main="",ylab="Ratio of means",xlab="",col="gray47",cex.lab=2,cex.axis=1.5)

title(main="(A) Scaffold 1631 (Z)",adj=0,cex.main=2.5)

abline(h=1,lwd=3,col="green3")

abline(h=2,lwd=3,col="orangered")

#plot2

par(mar=c(5,7,3,1))

plot(males_winmeans28/fem_winmeans28,type='l', main="",ylab="",xlab="",col="gray47",cex.lab=2,cex.axis=1.5)

title(main="(B) Scaffold 1628",adj=0,cex.main=2.5)

abline(h=1,lwd=3,col="green3")

abline(h=2,lwd=3,col="orangered")

#plot3

par(mar=c(5,7,3,1))

plot(males_winmeans45/fem_winmeans45,type='l', main="",ylab="Ratio of means", xlab="Distance (10kb)",col="gray47",cex.lab=2,cex.axis=1.5)

title(main="(C) Scaffold 1645",adj=0,cex.main=2.5)

abline(h=1,lwd=3,col="green3")

abline(h=2,lwd=3,col="orangered")

#plot4

par(mar=c(5,7,3,1))

plot(males_winmeans46/fem_winmeans46,type='l', main="",ylab="",xlab="Distance (10kb)",col="gray47",cex.lab=2,cex.axis=1.5)

title(main="(D) Scaffold 1646",adj=0,cex.main=2.5)

abline(h=1,lwd=3,col="green3")

abline(h=2,lwd=3,col="orangered")

dev.off()