Load libraries

library(kableExtra)
library(tidyverse)
library(edgeR)
library(WGCNA)
library(DESeq2)

Sample data

Supplementary Table M 1: Taiwan habu specimens

library individual day sex svl mass mix tissue
Pm_1_heart Pm_1 1 F 77.5 107.5 1 heart
Pm_10_heart Pm_10 2 F 63.4 82.2 1 heart
Pm_3_heart Pm_3 1 M 62.8 65.5 2 heart
Pm_5_heart Pm_5 1 M 76.5 136.9 1 heart
Pm_8_heart Pm_8 2 M 69.0 89.0 2 heart
Pm_1_kidney Pm_1 1 F 77.5 107.5 2 kidney
Pm_10_kidney Pm_10 2 F 63.4 82.2 2 kidney
Pm_3_kidney Pm_3 1 M 62.8 65.5 1 kidney
Pm_5_kidney Pm_5 1 M 76.5 136.9 2 kidney
Pm_8_kidney Pm_8 2 M 69.0 89.0 1 kidney
Pm_1_liver Pm_1 1 F 77.5 107.5 1 liver
Pm_10_liver Pm_10 2 F 63.4 82.2 1 liver
Pm_3_liver Pm_3 1 M 62.8 65.5 2 liver
Pm_5_liver Pm_5 1 M 76.5 136.9 1 liver
Pm_8_liver Pm_8 2 M 69.0 89.0 2 liver
Pm_1 Pm_1 1 F 77.5 107.5 1 venom gland
Pm_10 Pm_10 2 F 63.4 82.2 2 venom gland
Pm_11 Pm_11 2 F 76.0 122.1 1 venom gland
Pm_12 Pm_12 2 M 98.8 261.2 2 venom gland
Pm_13 Pm_13 2 F 61.5 54.6 2 venom gland
Pm_14 Pm_14 4 M 56.3 44.5 1 venom gland
Pm_15 Pm_15 4 F 69.5 86.7 2 venom gland
Pm_16 Pm_16 4 M 79.3 119.4 1 venom gland
Pm_17 Pm_17 4 F 68.7 73.9 2 venom gland
Pm_18 Pm_18 4 M 61.8 64.3 1 venom gland
Pm_19 Pm_19 4 F 74.0 91.6 1 venom gland
Pm_2 Pm_2 1 M 62.2 64.4 2 venom gland
Pm_20 Pm_20 4 F 79.9 128.9 2 venom gland
Pm_21 Pm_21 4 F 71.3 109.4 1 venom gland
Pm_22 Pm_22 8 F 68.2 79.9 2 venom gland
Pm_23 Pm_23 8 F 58.1 55.0 1 venom gland
Pm_24 Pm_24 8 F 69.0 101.9 1 venom gland
Pm_25 Pm_25 8 M 63.1 74.4 2 venom gland
Pm_26 Pm_26 8 F 65.4 74.1 1 venom gland
Pm_27 Pm_27 8 F 73.7 119.1 2 venom gland
Pm_28 Pm_28 8 M 69.6 88.0 2 venom gland
Pm_29 Pm_29 8 M 65.3 71.8 1 venom gland
Pm_3 Pm_3 1 M 62.8 65.5 2 venom gland
Pm_30 Pm_30 2 F 70.0 96.6 2 venom gland
Pm_4 Pm_4 1 F 77.0 97.6 1 venom gland
Pm_5 Pm_5 1 M 76.5 136.9 2 venom gland
Pm_6 Pm_6 1 F 78.7 149.5 1 venom gland
Pm_7 Pm_7 2 F 68.1 87.0 1 venom gland
Pm_8 Pm_8 2 M 69.0 89.0 1 venom gland
Pm_9 Pm_9 2 F 109.3 295.6 2 venom gland

Quality control

ERCC92 was added to the RNA-seq libraries for quality checking

fmt <- function() function(x) as.character(round(x,2))
rsem <- read_csv("./Data_sets/rsem.csv.gz",col_types = "ccnncnccc")
ercc <- read.csv("./Data_sets/ercc.csv.gz")

Supplementary Figure M 1: Spike in tpm

ggplot(ercc, aes(conc, tpm, color=library, group=library))+geom_point()+scale_x_continuous(labels = fmt(),trans="log2")+ scale_y_continuous(labels = fmt(),trans="log2")+theme_bw()+ xlab("Spike-in concentration (attamol/ul)")+ylab("Spike-in TPM") + geom_hline(yintercept = .05, color = "red")+guides(color=F)

There is generally good log-linear response between observed and expected values, though we can't detect anything lower than 0.05 TPM. Looks like one library got less spike in than the rest.

factors <- read_csv("./Data_sets/factors.csv", col_types = cols()) %>% filter(tissue == "venom gland")
rsemVenom <- rsem[grepl('Pm_\\d+$', rsem$library),]
# remove low expressed loci
keep <- rsemVenom %>% group_by(gene) %>% summarize(tpm = mean(tpm)) %>% filter(tpm>0.05) %>% select(gene) # genes to keep
## `summarise()` ungrouping output (override with `.groups` argument)
tpmVenom <- rsemVenom %>% select(library, gene, count, description, class) %>% filter(gene %in% keep$gene) %>%  spread(key = library, value = count)  %>%  remove_rownames %>% column_to_rownames("gene")
# do differential gene expression on venom components only
countsVenom <- rsemVenom %>% filter(!is.na(class)) %>% select(library, gene, count) %>%  spread(key = library, value = count) %>% 
  remove_rownames %>% column_to_rownames("gene")
#design matrix
design <- model.matrix(~factor(day)+factor(sex)+0, factors)
rownames(design) <- rownames(factors)
dge <- DGEList(countsVenom)
dge <- calcNormFactors(dge)
dat <- estimateGLMCommonDisp(dge, design)
dat <- estimateGLMTrendedDisp(dat, design)
dat <- estimateGLMTagwiseDisp(dat, design)
fit <- glmFit(dat,design)
mds <- plotMDS(dge)

mdsPlot <- data.frame(mds1=mds$cmdscale.out[,1], mds2=mds$cmdscale.out[,2], ID=rownames(mds$cmdscale.out),days=factors$day,sex=factors$sex,mass=factors$mass)

Supplementary Figure M 2: MDS mass and sex

ggplot(mdsPlot,aes(mds1,mds2,size=mass,shape=sex,color=as.factor(days)))+geom_point()+scale_colour_grey(start=0.8, end = 0)

maleVsfemale <- glmLRT(fit,coef=5)  
summary(caste.de <- decideTestsDGE(maleVsfemale, p=0.05, adjust="BH"))
##        factor(sex)M
## Down              0
## NotSig           76
## Up                0

Supplementary Figure M 3: Differences in venom expression between sexes is not significant

ggplot(na.omit(subset(rsem, tissue == "venom gland")), aes(class,tpm, group=sex)) + scale_y_log10()+ stat_summary(fun.data = "mean_cl_boot", position = position_dodge(width = .75), aes(colour = sex))

t.test(mass~sex,factors)
## 
##  Welch Two Sample t-test
## 
## data:  mass by sex
## t = 0.36085, df = 18.552, p-value = 0.7223
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -37.66358  53.32483
## sample estimates:
## mean in group F mean in group M 
##       105.95789        98.12727
t.test(svl~sex,factors)
## 
##  Welch Two Sample t-test
## 
## data:  svl by sex
## t = 0.71199, df = 19.303, p-value = 0.485
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -5.957906 12.111016
## sample estimates:
## mean in group F mean in group M 
##        72.59474        69.51818

Calculating soft-threshold

The lowest power to attain scale-free topology was 13 (R.sq = 0.903). Thereofre, we'll use a soft threshold of 13 to construct our network.

datExpr <- t(rlog(as.matrix(round(tpmVenom %>% select(-c(Pm_7, description,class))))))
## converting counts to integer mode
powers = c(seq(from = 10, to=16, by=1))

#call network topology analyis function
sft <- pickSoftThreshold(datExpr, powerVector = powers, verbose = 5, networkType = "signed")
## pickSoftThreshold: will use block size 2443.
##  pickSoftThreshold: calculating connectivity for given powers...
##    ..working on genes 1 through 2443 of 18313
##    ..working on genes 2444 through 4886 of 18313
##    ..working on genes 4887 through 7329 of 18313
##    ..working on genes 7330 through 9772 of 18313
##    ..working on genes 9773 through 12215 of 18313
##    ..working on genes 12216 through 14658 of 18313
##    ..working on genes 14659 through 17101 of 18313
##    ..working on genes 17102 through 18313 of 18313
##   Power SFT.R.sq slope truncated.R.sq mean.k. median.k. max.k.
## 1    10    0.865 -2.06          0.972   189.0     136.0    821
## 2    11    0.881 -2.01          0.978   143.0      97.0    727
## 3    12    0.897 -1.96          0.983   111.0      69.8    650
## 4    13    0.903 -1.92          0.981    87.3      50.9    585
## 5    14    0.911 -1.89          0.983    69.8      37.5    530
## 6    15    0.907 -1.86          0.976    56.6      27.9    483
## 7    16    0.913 -1.83          0.977    46.4      20.9    442
sizeGrWindow(9, 5)
par(mfrow = c(1,2))
cex1 = 0.9

Supplementary Figure M 4: Soft threshold to attain scale-free topology

# Scale-free topology fit index as a function of the soft-thresholding power
plot(sft$fitIndices[,1], -sign(sft$fitIndices[,3])*sft$fitIndices[,2],
     xlab="Soft Threshold (power)",ylab="Scale Free Topology Model Fit,signed R^2",type="n",
     main = paste("Scale independence"))
text(sft$fitIndices[,1], -sign(sft$fitIndices[,3])*sft$fitIndices[,2],
     labels=powers,cex=cex1,col="red")

# this line corresponds to using an R^2 cut-off of h
abline(h=0.90,col="red") 

plot(sft$fitIndices[,1], sft$fitIndices[,5],
     xlab="Soft Threshold (power)",ylab="Mean Connectivity", type="n",
     main = paste("Mean connectivity"))
text(sft$fitIndices[,1], sft$fitIndices[,5], labels=powers, cex=cex1,col="red")

Public RNA data

Supplementary Table M 2: Public RNA datasets

samples<-read.csv("./Data_sets/SRA_datasets.csv")
kable(cbind(samples)) %>%
  kable_styling() %>%
  scroll_box(width = "700px", height = "500px")
Taxa Tissue Library Project Layout Study X
Anole heart SRR579558 PRJNA176589 Paired Barbosa-Morais NL et al.,
Anole heart SRR5412152 PRJNA381064 Single Marin R et al.,
Anole heart SRR5412154 PRJNA381064 Single Marin R et al.,
Anole heart SRR5412155 PRJNA381064 Single Marin R et al.,
Anole liver SRR5412164 PRJNA381064 Single Marin R et al.,
Anole liver SRR5412166 PRJNA381064 Single Marin R et al.,
Anole liver SRR5412167 PRJNA381064 Single Marin R et al.,
Anole kidney SRR579557 PRJNA176589 Paired Barbosa-Morais NL et al.,
Anole kidney SRR5412157 PRJNA381064 Single Marin R et al.,
Anole kidney SRR5412158 PRJNA381064 Single Marin R et al.,
Anole kidney SRR5412161 PRJNA381064 Single Marin R et al.,
Canis heart SRR8996963 PRJNA396033 Paired https://www.broadinstitute.org/mammals-models/broad-institute-canine-genomic-resources
Canis heart SRR6206901 PRJNA396033 Paired https://www.broadinstitute.org/mammals-models/broad-institute-canine-genomic-resources
Canis heart SRR8996973 PRJNA396033 Paired https://www.broadinstitute.org/mammals-models/broad-institute-canine-genomic-resources
Canis liver SRR8997009 PRJNA396033 Paired https://www.broadinstitute.org/mammals-models/broad-institute-canine-genomic-resources
Canis liver SRR8997042 PRJNA396033 Paired https://www.broadinstitute.org/mammals-models/broad-institute-canine-genomic-resources
Canis liver SRR6206911 PRJNA415552 Paired Chen J et al.
Canis kidney SRR8996953 PRJNA396033 Paired https://www.broadinstitute.org/mammals-models/broad-institute-canine-genomic-resources
Canis kidney SRR8997002 PRJNA396033 Paired https://www.broadinstitute.org/mammals-models/broad-institute-canine-genomic-resources
Canis kidney SRR6206906 PRJNA396033 Paired https://www.broadinstitute.org/mammals-models/broad-institute-canine-genomic-resources
Canis salivary gland SRR5889344 PRJNA396033 Paired https://www.broadinstitute.org/mammals-models/broad-institute-canine-genomic-resources
Canis salivary gland SRR8997053 PRJNA396033 Paired https://www.broadinstitute.org/mammals-models/broad-institute-canine-genomic-resources
Gal heart SRR579554 PRJNA176589 Paired Barbosa-Morais NL et al.,
Gal heart SRR5412246 PRJNA381064 Single Marin R et al.,
Gal heart SRR5412248 PRJNA381064 Single Marin R et al.,
Gal liver SRR579552 PRJNA176589 Paired Barbosa-Morais NL et al.,
Gal liver SRR5412254 PRJNA381064 Single Marin R et al.,
Gal liver SRR5412256 PRJNA381064 Single Marin R et al.,
Gal kidney SRR5412256 PRJNA381064 Single Marin R et al.,
Gal kidney SRR5412250 PRJNA381064 Single Marin R et al.,
Gal kidney SRR5412251 PRJNA381064 Single Marin R et al.,
Homo heart SRR306847 PRJNA143627 Single Brawand D et al.,
Homo heart SRR306848 PRJNA143627 Single Brawand D et al.,
Homo heart SRR306849 PRJNA143627 Single Brawand D et al.,
Homo liver SRR306856 PRJNA143627 Single Brawand D et al.,
Homo liver SRR5412182 PRJNA381064 Single Marin R et al.,
Homo liver SRR5412183 PRJNA381064 Single Marin R et al.,
Homo kidney SRR306851 PRJNA143627 Single Brawand D et al.,
Homo kidney SRR306852 PRJNA143627 Single Brawand D et al.,
Homo kidney SRR5412179 PRJNA381064 Single Marin R et al.,
Homo salivary gland SRR10246417 PRJNA576455 Paired Alame M et al.
Homo salivary gland SRR10246418 PRJNA576455 Paired Alame M et al.
Mus heart SRR453173 SRP012040 Paired GSE36025: Long RNA-seq from ENCODE/Cold Spring Harbor Lab
Mus heart SRR453174 SRP012040 Paired GSE36025: Long RNA-seq from ENCODE/Cold Spring Harbor Lab
Mus heart SRR453175 SRP012040 Paired GSE36025: Long RNA-seq from ENCODE/Cold Spring Harbor Lab
Mus liver SRR453152 SRP012040 Paired GSE36025: Long RNA-seq from ENCODE/Cold Spring Harbor Lab
Mus liver SRR453153 SRP012040 Paired GSE36025: Long RNA-seq from ENCODE/Cold Spring Harbor Lab
Mus liver SRR453154 SRP012040 Paired GSE36025: Long RNA-seq from ENCODE/Cold Spring Harbor Lab
Mus liver SRR453155 SRP012040 Paired GSE36025: Long RNA-seq from ENCODE/Cold Spring Harbor Lab
Mus kidney SRR453147 SRP012040 Paired GSE36025: Long RNA-seq from ENCODE/Cold Spring Harbor Lab
Mus kidney SRR453148 SRP012040 Paired GSE36025: Long RNA-seq from ENCODE/Cold Spring Harbor Lab
Mus kidney SRR453149 SRP012040 Paired GSE36025: Long RNA-seq from ENCODE/Cold Spring Harbor Lab
Mus salivary gland SRR8133129 PRJNA501787 Paired https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE121954 https://www.ncbi.nlm.nih.gov/sra/SRX4954232[accn]
Mus salivary gland SRR8133134 PRJNA501787 Paired https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE121954 https://www.ncbi.nlm.nih.gov/sra/SRX4954232[accn]
Mus salivary gland SRR8133136 PRJNA501787 Paired https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE121954 https://www.ncbi.nlm.nih.gov/sra/SRX4954232[accn]
Mus salivary gland SRR8133137 PRJNA501787 Paired https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE121954 https://www.ncbi.nlm.nih.gov/sra/SRX4954232[accn]
Naja liver SRR8754967 PRJNA527614 Single Suryamohan K
Naja liver SRR8754989 PRJNA527614 Single Suryamohan K
Naja kidney SRR8754966 PRJNA527614 Single Suryamohan K
Naja kidney SRR8754970 PRJNA527614 Single Suryamohan K
Naja salivary gland SRR8754972 PRJNA527614 Single Suryamohan K
Naja venom gland SRR8754964 PRJNA527614 Single Suryamohan K
Naja venom gland SRR8754977 PRJNA527614 Single Suryamohan K
Naja venom gland SRR8754980 PRJNA527614 Single Suryamohan K
Naja venom gland SRR8754981 PRJNA527614 Single Suryamohan K
Naja venom gland SRR8754986 PRJNA527614 Single Suryamohan K
Naja venom gland SRR8754988 PRJNA527614 Single Suryamohan K
Pan heart SRR306819 PRJNA143627 Single Brawand D et al.,
Pan heart SRR306820 PRJNA143627 Single Brawand D et al.,
Pan liver SRR306823 PRJNA143627 Single Brawand D et al.,
Pan kidney SRR306821 PRJNA143627 Single Brawand D et al.,
Pan kidney SRR306822 PRJNA143627 Single Brawand D et al.,
Xeno heart SRR579562 PRJNA176589 Paired Barbosa-Morais NL et al.,
Xeno heart SRR5412268 PRJNA381064 Single Marin R et al.,
Xeno liver SRR579561 PRJNA176589 Paired Barbosa-Morais NL et al.,
Xeno liver SRR5412274 PRJNA381064 Single Marin R et al.,
Xeno liver SRR5412276 PRJNA381064 Single Marin R et al.,
Xeno kidney SRR579562 PRJNA176589 Paired Barbosa-Morais NL et al.,
Xeno kidney SRR5412270 PRJNA381064 Single Marin R et al.,
Xeno kidney SRR5412271 PRJNA381064 Single Marin R et al.,

CAFE error model and output

Supplementary Table M 3: CAFE error model

cat ./CAFE/res/Base_error_model.txt
maxcnt: 1
cntdiff: -1 0 1
0 0 0.831338 0.168662
1 0.168662 0.662677 0.168662

Two lambda CAFE output (lambda1 and lambda2)

cat ./CAFE/res/two_lambdamodel_results.txt
Model Base Result: 6909.73
Lambda: 0.0017969444427155, 0.0064501608736268
Epsilon: 0.168662

Supplementary Figure M 5: Null distribution of likelihood rations obtained from simulation in CAFE

Once the log-likelihoods were obtained, lhtest.R script (provided by CAFE) was used to create histogram with a null distribution obtained from simulations. Significance is determined by how far left the observed likelihood ratio (2 x lnLglobal - lnLmulti) would fall on the tail of the distribution. In our case the likelihood ratio count would fall on the far left of the distribution indicating a very low p value