Post date: Aug 23, 2019 10:57:22 PM
I basically repeated the mapping analyses with gemma described here (T. knulli, mapping performance on RW versus C: initial results and PCA by LG), but used weight change between 15 and 21 days as the phenotype.
Here is the new code in formatPhenoGeno.R for this:
## weight change
C<-which(ph$Treatment.March.19=="C")
RW<-which(ph$Treatment.March.19=="RW")
phC<-ph[C,]
phRW<-ph[RW,]
gemma_phC2<-matrix(NA,nrow=length(C),ncol=2)
lmo<-lm(phC$Weight.at.day.21-phC$Weight.at.Day.15 ~ phC$Stage.at.Day.15 * phC$Sex.at.Day.21)
gemma_phC2[is.na(phC$Weight.at.day.21)==FALSE,1]<-lmo$residuals
lmo<-lm(phC$Weight.at.day.21-phC$Weight.at.Day.15 ~ phC$Sex.at.Day.21)
gemma_phC2[is.na(phC$Weight.at.day.21)==FALSE,2]<-lmo$residuals
gemma_phRW2<-matrix(NA,nrow=length(RW),ncol=2)
lmo<-lm(phRW$Weight.at.day.21-phRW$Weight.at.Day.15 ~ phRW$Stage.at.Day.15 * phRW$Sex.at.Day.21)
gemma_phRW2[is.na(phRW$Weight.at.day.21)==FALSE,1]<-lmo$residuals
lmo<-lm(phRW$Weight.at.day.21-phRW$Weight.at.Day.15 ~ phRW$Sex.at.Day.21)
gemma_phRW2[is.na(phRW$Weight.at.day.21)==FALSE,2]<-lmo$residuals
write.table(round(gemma_phC2,5),file="pheno_C_dw.txt",row.names=FALSE,col.names=FALSE,quote=FALSE)
write.table(round(gemma_phRW2,5),file="pheno_RW_dw.txt",row.names=FALSE,col.names=FALSE,quote=FALSE)
save(list=ls(),file="gp.rdat")
Note that phenotype 1 controls for stage (at day 15) whereas phenotype 2 does not. I then ran the gemma BSLMM as before:
#!/usr/bin/perl
#
# run gemma jobs
#
use Parallel::ForkManager;
my $max = 20;
my $pm = Parallel::ForkManager->new($max);
foreach $g (@ARGV){
$ph = $g;
$ph =~ s/mod_g/ph/;
$ph =~ s/\.txt/_dw.txt/;
$base = $g;
$base =~ s/mod_geno_//;
$base =~ s/\.txt//;
foreach $ch (0..9){
system "sleep 2\n";
foreach $trait (1..2){
#foreach $trait (1..3){
$pm->start and next; ## fork
$out = "o_$base"."_dw_ph$trait"."_ch$ch";
#$out = "o_$base"."_ph$trait"."_ch$ch";
system "gemma -g $g -p $ph -bslmm 1 -n $trait -maf 0.0 -w 200000 -s 1000000 -o $out\n";
$pm->finish;
}
}
}
$pm->wait_all_children;
The results are in output_bslmm_change. I plotted summaries with plots.R (see below). I don't really see much of a signal at all. Perhaps we stick with weight.
files<-list.files(pattern="mav")
L<-4
dat<-vector("list",L)
for(i in 1:L){
dat[[i]]<-scan(files[i])
}
snps<-read.table("snpList.txt",header=F)
csn<-as.numeric(snps[,1]==12033)
tits<-c("Weight change 1, C","Weight change 2, C",
"Weight change 1, RW","Weight change 2, RW")
pdf("TknulBSLMM_mav_dw.pdf",width=8,height=8)
par(mfrow=c(2,1))
par(mar=c(4,5,2,1))
for(i in 1:L){
plot(abs(dat[[i]]),pch=19,col=c("darkgray","cadetblue")[csn+1],xlab="SNP",ylab="Mod. average eff.",cex.lab=1.4)
title(main=tits[i])
}
dev.off()
##########
files<-list.files(pattern="pip_")
L<-4
dat<-vector("list",L)
for(i in 1:L){
dat[[i]]<-scan(files[i])
}
snps<-read.table("snpList.txt",header=F)
tits<-c("Weight change 1, C","Weight change 2, C",
"Weight change 1, RW","Weight change 2, RW")
pdf("TknulBSLMM_QTL_dw.pdf",width=9,height=9)
par(mfrow=c(2,2))
par(mar=c(4,5,2,1))
n<-table(snps[,1])
for(i in 1:L){
qtl<-tapply(X=dat[[i]],INDEX=snps[,1],sum)
cs<-rep("darkgray",length(qtl))
cs[which(names(qtl)==12033)]<-"cadetblue"
plot(as.numeric(n),qtl,pch=19,col=cs,xlab="SNP",ylab="PIP",cex.lab=1.4)
title(main=tits[i])
}
dev.off()
##########
files<-list.files(pattern="pip_")
L<-4
dat<-vector("list",L)
for(i in 1:L){
dat[[i]]<-scan(files[i])
}
snps<-read.table("snpList.txt",header=F)
csn<-as.numeric(snps[,1]==12033)
tits<-c("Weight change 1, C","Weight change 2, C",
"Weight change 1, RW","Weight change 2, RW")
pdf("TknulBSLMM_dw_pip.pdf",width=8,height=8)
par(mfrow=c(2,1))
par(mar=c(4,5,2,1))
for(i in 1:L){
plot(dat[[i]],pch=19,col=c("darkgray","cadetblue")[csn+1],xlab="SNP",ylab="PIP.",cex.lab=1.4)
title(main=tits[i])
}
dev.off()