SLICK Holstein mammary gland RNA-Seq methods

Bos taurus reference

Genome assembly: ARS-UCD1.3

Download using the NCBI datasets tool

datasets download genome accession GCF_002263795.2 \
    --include gff3,rna,cds,protein,genome,seq-report \
    --filename $DIR/0-ref/GCF_002263795.2.zip

extract

cd $DIR/0-ref
unzip GCF_002263795.2.zip

create salmon database

PREFIX=$DIR/0-ref/ncbi_dataset/data/GCF_002263795.2

cat $PREFIX/rna.fna $PREFIX/GCF_002263795.2_ARS-UCD1.3_genomic.fna \
    > $DIR/0-ref/bovine.db.fa
sed -n '/>/s/>\([^ ]*\).*/\1/p' $PREFIX/GCF_002263795.2_ARS-UCD1.3_genomic.fna \
    > $DIR/0-ref/decoys.txt

$DIR/apps/salmon-1.9.0_linux_x86_64/bin/salmon index \
    -t $DIR/0-ref/bovine.db.fa \
    -i $DIR/0-ref/bovine.db \
    --decoys $DIR/0-ref/decoys.txt \
    -p 48

Raw data

Library	Forward	Reverse
SLICK_2041	SLICK_2041_1.fq.gz	SLICK_2041_2.fq.gz
SLICK_2128	SLICK_2128_1.fq.gz	SLICK_2128_2.fq.gz
SLICK_2129	SLICK_2129_1.fq.gz	SLICK_2129_2.fq.gz
SLICK_2171	SLICK_2171_1.fq.gz	SLICK_2171_2.fq.gz
SLICK_2202	SLICK_2202_1.fq.gz	SLICK_2202_2.fq.gz
SLICK_2287	SLICK_2287_1.fq.gz	SLICK_2287_2.fq.gz
SLICK_466	SLICK_466_1.fq.gz	SLICK_466_2.fq.gz
WT_2058	WT_2058_1.fq.gz	WT_2058_2.fq.gz
WT_2068	WT_2068_1.fq.gz	WT_2068_2.fq.gz
WT_2089	WT_2089_1.fq.gz	WT_2089_2.fq.gz
WT_2116	WT_2116_1.fq.gz	WT_2116_2.fq.gz
WT_2179	WT_2179_1.fq.gz	WT_2179_2.fq.gz
WT_2226	WT_2226_1.fq.gz	WT_2226_2.fq.gz
WT_2273	WT_2273_1.fq.gz	WT_2273_2.fq.gz

Salmon

RAW=$DIR/1-raw/

for name in "${!data[@]}"; do
  readarray -t lib <<<"${data[$name]}"

  $DIR/apps/salmon-1.9.0_linux_x86_64/bin/salmon quant \
      -i $DIR/0-ref/bovine.db \
      -l A \
      -1 $RAW/${lib[0]} \
      -2 $RAW/${lib[1]} \
      --validateMappings \
      -o $DIR/2-salmon/$name \
      --threads 48
done

awk '/^Observed/ {obs=$2}
     /Mapping rate/ {rate=$NF}
     FNR == 1 {
         if(NR!=FNR) printf "%s\t%'"'"'d\t%0.2f%%\n", name, obs, rate;
         obs_all += obs;
         rate_all += rate;
         c++;
         name=FILENAME;
         sub(".*/2-salmon/", "", name);
         sub("/.*", "", name);
          }
     END {
         obs_all += obs;
         rate_all += rate;
         c++;

         printf "%s\t%'"'"'d\t%0.2f%%\n", name, obs, rate;
         printf "%s\t%'"'"'d\t%0.2f%%\n", "Average", obs_all/c, rate_all/c;
         }' OFS="\t" \
         $(printf "$DIR/2-salmon/%s/logs/salmon_quant.log " "${names[@]}")

Library	Observed Fragments	Mapping Rate
SLICK_2041	30,023,454	84.85%
SLICK_2128	34,213,368	89.53%
SLICK_2129	30,141,707	90.01%
SLICK_2171	32,747,566	81.12%
SLICK_2202	39,783,888	88.73%
SLICK_2287	35,395,347	90.13%
SLICK_466	31,035,608	87.05%
WT_2058	34,627,012	86.17%
WT_2068	31,644,602	88.24%
WT_2089	35,163,794	83.39%
WT_2116	35,970,980	89.30%
WT_2179	30,755,764	90.12%
WT_2226	40,432,327	89.07%
WT_2273	29,421,735	89.04%
Average	31,423,810	87.62%

EdgeR

Load data

library(edgeR)
library(statmod)
library(tximport)
library(data.table)
library(tidyverse)
library(ggrepel)

meta.data <- data.frame(library=c("SLICK_2041", "SLICK_2128", "SLICK_2129",
                                  "SLICK_2171", "SLICK_2202", "SLICK_2287",
                                  "SLICK_466", "WT_2058", "WT_2068", "WT_2089",
                                  "WT_2116", "WT_2179", "WT_2226", "WT_2273"),
                        group=factor(rep(c("SL", "WT"), each=7), levels=c("WT", "SL")))

files <- file.path('2-salmon', meta.data$library,'quant.sf')
names(files) <- meta.data$library


design <- model.matrix(~ group, data=meta.data)
colnames(design) <- make.names(gsub("group", "", colnames(design)))

## Create tx2gene from gff
txdb <- fread("0-ref/ncbi_dataset/data/GCF_002263795.2/genomic.gff",
              skip=8, fill=T, header=F, sep="\t")[V3=='mRNA',] %>%
  dplyr::select('Attr'=V9) %>%
  as.data.frame %>%
  tibble::rowid_to_column(var='rowid') %>%
  separate_rows(Attr, sep=';') %>%
  separate(Attr, into=c('key', 'value'), sep='=') %>%
  spread(key, value) %>%
  mutate(GeneID = sub(".*GeneID:([^,]*).*", "\\1", Dbxref))

tx2gene <- dplyr::select(txdb, Name, GeneID)

txi <- tximport(files, type = "salmon", tx2gene = tx2gene)

Normalize counts, Filter Genes, and Estimate Dispersions

cts <- txi$counts
normMat <- txi$length

## Obtaining per-observation scaling factors for length, adjusted to avoid
## changing the magnitude of the counts.
normMat <- normMat/exp(rowMeans(log(normMat)))
normCts <- cts/normMat

 # Computing effective library sizes from scaled counts, to account for
# composition biases between samples.
eff.lib <- calcNormFactors(normCts) * colSums(normCts)

# Combining effective library sizes with the length factors, and calculating
# offsets for a log-link GLM.
normMat <- sweep(normMat, 2, eff.lib, "*")
normMat <- log(normMat)

# Creating a DGEList object for use in edgeR.
d <- DGEList(cts) %>% scaleOffset(normMat)

# filtering using the design information
keep <- filterByExpr(d, design)

y <- d[keep, ,keep.lib.sizes=F ] %>%
  calcNormFactors %>%
  estimateDisp(design, robust=T)

unique genes

as.data.frame(txi$counts) %>%
  rownames_to_column('GeneID') %>%
gather(-GeneID, key="library", value="count") %>%
  left_join(meta.data) %>%
  filter(count > 0) %>%
  group_by(GeneID, group) %>%
  count %>%
  spread(group, n) %>%
  filter(is.na(WT) | is.na(SL)) %>%
  gather(-GeneID, key="Type", value="Samples") %>%
  filter(Samples >= 4)

GeneID	Type	Samples
112441545	WT	6
112447300	WT	5
505979	WT	4
505993	WT	4
523830	WT	4
535831	WT	4
615310	WT	4
786611	WT	5
101905578	SL	4
101906363	SL	4
104968609	SL	5
112449108	SL	4
512850	SL	4
516736	SL	5
524709	SL	4
526745	SL	4
540297	SL	4
613925	SL	5
783189	SL	4
790886	SL	4

Unique Gene Venn Diagram

unique.gene.venn.plot <- as.data.frame(txi$counts) %>%
  rownames_to_column('gene') %>%
gather(-gene, key="library", value="count") %>%
  left_join(meta.data) %>%
  filter(count > 0) %>%
  group_by(gene, group) %>%
  count %>%
  spread(group, n) %>%
  filter(is.na(WT) | is.na(SL)) %>%
  gather(-gene, key="type", value="count") %>%
  filter(!is.na(count)) %>%
  arrange(count) %>%
  mutate(type=fct_rev(type)) %>%
  ggplot(aes(type,1,fill=count)) +
  geom_col() +
  scale_x_discrete(breaks=c("WT", "SL"), labels=c("Wild-Type", "SLICK")) +
  scale_y_continuous(expand=c(0,0)) +
  scale_fill_distiller(palette = "PuRd") +
  guides(fill=guide_colorbar(title="Number of \nSamples")) +
  labs(y="Number of unique genes", x="Hair Coat Type") +
  theme_bw()+
  theme(
    panel.grid = element_blank(),
    axis.text = element_text(size=14),,
    axis.title = element_text(size=16),
    legend.text = element_text(size=14),
    legend.title = element_text(size=16)
    )

unique.gene.venn.plot

Filtered/Surviving gene summary
```
table(keep)
    
```

MDS

mds.data <- plotMDS(y, plot=F)

mds.plot <- cbind(meta.data, x=mds.data$x, y=mds.data$y) %>%
  mutate(library = sub('.*_', '', library)) %>%
  ggplot(aes(x,y, color=group, label=library)) +
  geom_point() +
  geom_text_repel(show.legend = F) +
  scale_color_manual(values=c( "#cc546e", "#4db598"), labels=c("Wild-Type", "SLICK")) +
  xlab(sprintf("Dimension 1 (%0.2f%%)", mds.data$var.explained[1]*100)) +
  ylab(sprintf("Dimension 2 (%0.2f%%)", mds.data$var.explained[2]*100)) +
  labs(color="Hair Coat Type") +
  theme_bw() +
  theme(legend.position = c(0.8,0.8),
        legend.background = element_rect(fill='white', color="grey"),
        panel.grid = element_blank()) +
    theme( legend.position = c(0.8,0.9),
    panel.grid.minor = element_blank(),
    axis.text = element_text(size=14),,
    axis.title = element_text(size=16),
    legend.text = element_text(size=14),
    legend.title = element_text(size=16)
    )
mds.plot

Panel MDS/Venn

library(cowplot)

plot_grid(mds.plot, unique.gene.venn.plot, labels=c("A", "B"), rel_widths = c(3,2))

Mean-Difference (MD) plots

par(mfrow = c(3,5))
for(x in 1:14){
  plotMD(y, column=x, ylim=c(-8,8))
  abline(h=0, col="red", lty=2, lwd=2)
}

Biological Coefficient of Variation (BCV)
```
plotBCV(y)
    
```

DEG LRT test

library(ascii)
options(asciiType = "org", width = 80)

fit <- glmFit(y,design,robust=T)
lrt <- glmLRT(fit)

desc <- txdb %>%
  dplyr::select(GeneID, gene, product) %>%
  distinct(GeneID, .keep_all=T)

topTags(lrt, n=Inf) %>%
  as.data.frame %>%
  filter(PValue <= 0.01 & abs(logFC) > 0.6) %>%
  tibble::rownames_to_column("GeneID") %>%
  left_join(desc, by="GeneID")  %>%
  ascii(digits=4)

	GeneID	logFC	logCPM	LR	PValue	FDR	gene	product
1	513329	2.2834	1.2464	36.0838	0.0000	0.0000	LOC513329	major allergen Equ c 1%2C transcript variant X2
2	100139092	1.3766	0.3916	16.2967	0.0001	0.2886	GPR34	G protein-coupled receptor 34%2C transcript variant X2
3	511941	0.9035	3.1355	15.9611	0.0001	0.2886	CEACAM19	carcinoembryonic antigen related cell adhesion molecule 19
4	519145	1.5596	-0.5616	13.8251	0.0002	0.6717	LOC519145	beta-galactosidase-1-like protein 2%2C transcript variant X1
5	100336825	0.9564	0.9602	11.5818	0.0007	0.9998	TTLL10	tubulin tyrosine ligase like 10%2C transcript variant X2
6	280839	1.3017	0.4869	10.9832	0.0009	0.9998	LHB	luteinizing hormone beta polypeptide%2C transcript variant X1
7	100298265	1.6411	-0.0646	10.5104	0.0012	0.9998	CPAMD8	C3 and PZP like%2C alpha-2-macroglobulin domain containing 8%2C transcript variant X1
8	783257	-0.6126	2.5763	10.4299	0.0012	0.9998	ABCB6	ATP binding cassette subfamily B member 6
9	788425	1.3738	-0.2474	10.3593	0.0013	0.9998	LOC788425	aflatoxin B1 aldehyde reductase member 4%2C transcript variant X2
10	505551	1.2709	-1.0031	10.2327	0.0014	0.9998	MYCBPAP	MYCBP associated protein
11	539467	0.9888	3.1842	10.0680	0.0015	0.9998	PKDCC	protein kinase domain containing%2C cytoplasmic
12	286858	1.2357	8.3789	9.7174	0.0018	0.9998	PTGDS	prostaglandin D2 synthase%2C transcript variant X2
13	530346	1.2841	-0.0605	9.5219	0.0020	0.9998	ZNF621	zinc finger protein 621
14	519940	-1.5876	1.5054	8.7812	0.0030	0.9998	RDH16	retinol dehydrogenase 16 (all-trans)%2C transcript variant X5
15	100336551	-0.6278	2.7815	8.7279	0.0031	0.9998	ABCC9	ATP binding cassette subfamily C member 9%2C transcript variant X3
16	517459	-2.9386	2.1089	8.6482	0.0033	0.9998	GABRP	gamma-aminobutyric acid type A receptor pi subunit%2C transcript variant X2
17	286791	-1.9281	1.7573	8.5654	0.0034	0.9998	BDA20	major allergen BDA20%2C transcript variant X1
18	784038	-5.4453	5.0321	8.5339	0.0035	0.9998	LOC784038	testis-specific H1 histone
19	518047	-1.9013	-0.2802	8.4813	0.0036	0.9998	SLCO1C1	solute carrier organic anion transporter family member 1C1%2C transcript variant X2
20	534578	0.7251	3.5213	8.4583	0.0036	0.9998	LOC534578	vascular cell adhesion molecule 1-like
21	319095	-0.8635	4.0900	8.3819	0.0038	0.9998	ADCYAP1R1	ADCYAP receptor type I%2C transcript variant X5
22	524115	1.0508	-0.3196	7.9895	0.0047	0.9998	GABRE	gamma-aminobutyric acid type A receptor epsilon subunit
23	534512	1.1463	1.0113	7.8659	0.0050	0.9998	HMCN2	hemicentin 2%2C transcript variant X1
24	404129	1.4654	-0.6432	7.7895	0.0053	0.9998	DGAT2	diacylglycerol O-acyltransferase 2
25	100295741	1.3779	4.1745	7.7004	0.0055	0.9998	ZG16B	zymogen granule protein 16B
26	767895	0.6997	2.4466	7.6528	0.0057	0.9998	WDR60	WD repeat domain 60%2C transcript variant X4
27	534505	-1.8257	1.7249	7.6519	0.0057	0.9998	CPXM2	carboxypeptidase X%2C M14 family member 2
28	785756	-2.5650	3.6500	7.5531	0.0060	0.9998	LOC785756	androgen binding protein beta-like
29	539733	0.7298	2.1718	7.4913	0.0062	0.9998	RWDD2A	RWD domain containing 2A%2C transcript variant X1
30	613448	1.3919	0.7843	7.4564	0.0063	0.9998	AK5	adenylate kinase 5
31	515025	0.6691	3.0265	7.2853	0.0070	0.9998	BIN2	bridging integrator 2
32	523124	0.7462	0.7058	7.2466	0.0071	0.9998	PSD	pleckstrin and Sec7 domain containing%2C transcript variant X2
33	783399	4.5981	4.0281	6.8908	0.0087	0.9998	LOC783399	major allergen Equ c 1-like%2C transcript variant X1
34	615111	-0.6952	4.0500	6.8817	0.0087	0.9998	RNF150	ring finger protein 150%2C transcript variant X5
35	615896	-2.1402	1.3379	6.8148	0.0090	0.9998	SYCP3	synaptonemal complex protein 3%2C transcript variant X2
36	515653	0.6979	2.0786	6.7194	0.0095	0.9998	SGSM1	small G protein signaling modulator 1%2C transcript variant X1
37	527572	1.4034	0.7599	6.6427	0.0100	0.9998	GPR37L1	G protein-coupled receptor 37 like 1

MD plot

names <- dplyr::select(txdb, GeneID, gene) %>%
  distinct(GeneID, .keep_all=T)

topTags(lrt, n=Inf)$table %>%
    as.data.frame %>%
    rownames_to_column("GeneID") %>%
    left_join(names) %>%
    mutate(Sig = (PValue <= 0.01 & abs(logFC) >= 0.6),
          Dir = ifelse(Sig, ifelse(logFC < 0, "Down", "Up"), "Not Significant"),
          label=if_else(Sig, gene, NA)) %>%
     arrange(Sig) %>%
ggplot(aes(logCPM, logFC, color=Dir, size=Sig, label=label)) +
  geom_point() +
  geom_label_repel(size = 4, max.overlaps=100) +
  scale_size_manual(values=c(2/5, 1), guide=F) +
  scale_color_manual(values=c(Up="red", Down="blue", "grey"),
                     labels=c(Up="Up (25)", Down="Down (12)")) +
  xlab(expression("log"[2]*"CPM")) +
  ylab(expression("log"[2]*"FC")) +
  labs(color="Significant\n(p-value <= 0.01, |logFC| >= 0.6)") +
  theme_bw() +
  theme( legend.position = c(0.8,0.9),
    panel.grid.minor = element_blank(),
    axis.text = element_text(size=14),,
    axis.title = element_text(size=16),
    legend.text = element_text(size=14),
    legend.title = element_text(size=16)
    )

Volcano

topTags(lrt, n=Inf)$table %>%
   mutate(Sig = (PValue <= 0.01 & abs(logFC) >= 0.6),
          Dir = ifelse(Sig, ifelse(logFC < 0, "Down", "Up"), "Not Significant")) %>%
     arrange(Sig) %>%
ggplot(aes(logFC, -log(PValue,10), color=Dir, size=Sig)) + # -log10 conversion
  geom_point() +
  scale_size_manual(values=c(2/5, 1), guide=F) +
  scale_color_manual(values=c(Up="red", Down="blue", "grey")) +
  xlab(expression("log"[2]*"FC")) +
  ylab(expression("-log"[10]*"p-value")) +
  theme_minimal() +
  theme(legend.position = 'none')

KEGG

Download pathways and run fry

kg.list <- getGeneKEGGLinks("bta")
kg.sets <- split(kg.list$GeneID, kg.list$PathwayID)

kg.names <- getKEGGPathwayNames('bta')

kegg.fry <- fry(y, kg.sets, design)

Get pathways with p-value <= 0.05

kegg.fry %>%
  rownames_to_column('PathwayID') %>%
  mutate(PathwayID=sub('path:', '', PathwayID)) %>%
  left_join(kg.names, by="PathwayID") %>%
  filter(PValue <= 0.05) %>%
  mutate(Description = sub(' - Bos taurus \\(cow\\)', '', Description)) %>%
  ascii()

	PathwayID	NGenes	Direction	PValue	FDR	PValue.Mixed	FDR.Mixed	Description
1	bta00515	21.00	Up	0.01	1.00	0.42	0.99	Mannose type O-glycan biosynthesis
2	bta00590	42.00	Up	0.01	1.00	0.28	0.99	Arachidonic acid metabolism
3	bta00591	16.00	Up	0.04	1.00	0.27	0.99	Linoleic acid metabolism
4	bta00340	15.00	Down	0.04	1.00	0.31	0.99	Histidine metabolism

bta00515

Mannose type O-glycan biosynthesis

library(pathview)
pathway="00590"
pathview(pathway.id=pathway,
         kegg.dir = '4-kegg/',
         gene.data = setNames(lrt$table$logFC, row.names(lrt$table)),
         species='bta')

bta00590

Arachidonic acid metabolism

bta00591

Linoleic acid metabolism

bta00340

Histidine metabolism

GO

library(org.Bt.eg.db)
library(GO.db)

go.db <- select(org.Bt.eg.db,
       keys=keys(org.Bt.eg.db, keytype = "ENTREZID"),
       columns="GO",
       keytype="ENTREZID") %>%
  filter(!is.na(GO))


go.sets <- split(go.db$ENTREZID, go.db$GO)

go.names <- select(GO.db,
                   keys=names(go.sets),
                   columns=c("GOID", "TERM", "ONTOLOGY"),
                   keytype="GOID")

go.fry <- fry(y, go.sets, design)

Get go terms with p-value <= 0.05

go.fry %>%
  rownames_to_column('GOID') %>%
  filter(PValue <= 0.05) %>%
  left_join(go.names, by='GOID') %>%
  dplyr::select(-ends_with('.Mixed')) %>%
  ascii()

	GOID	NGenes	Direction	PValue	FDR	TERM	ONTOLOGY
1	GO:1903598	3.00	Down	0.00	1.00	positive regulation of gap junction assembly	BP
2	GO:0046628	1.00	Down	0.00	1.00	positive regulation of insulin receptor signaling pathway	BP
3	GO:0051302	6.00	Up	0.00	1.00	regulation of cell division	BP
4	GO:0005643	5.00	Down	0.00	1.00	nuclear pore	CC
5	GO:0005388	2.00	Down	0.00	1.00	P-type calcium transporter activity	MF
6	GO:0008265	1.00	Down	0.00	1.00	Mo-molybdopterin cofactor sulfurase activity	MF
7	GO:0043545	1.00	Down	0.00	1.00	molybdopterin cofactor metabolic process	BP
8	GO:0102867	1.00	Down	0.00	1.00	molybdenum cofactor sulfurtransferase activity	MF
9	GO:0140768	2.00	Down	0.00	1.00	protein ADP-ribosyltransferase-substrate adaptor activity	MF
10	GO:0050693	2.00	Down	0.00	1.00	LBD domain binding	MF
11	GO:0030151	2.00	Down	0.00	1.00	molybdenum ion binding	MF
12	GO:0021740	1.00	Down	0.00	1.00	principal sensory nucleus of trigeminal nerve development	BP
13	GO:0021960	1.00	Down	0.00	1.00	anterior commissure morphogenesis	BP
14	GO:0044300	1.00	Down	0.00	1.00	cerebellar mossy fiber	CC
15	GO:0060486	1.00	Down	0.00	1.00	club cell differentiation	BP
16	GO:0060509	1.00	Down	0.00	1.00	type I pneumocyte differentiation	BP
17	GO:0061141	1.00	Down	0.00	1.00	lung ciliated cell differentiation	BP
18	GO:0071679	1.00	Down	0.00	1.00	commissural neuron axon guidance	BP
19	GO:2000791	1.00	Down	0.00	1.00	negative regulation of mesenchymal cell proliferation involved in lung development	BP
20	GO:2000795	1.00	Down	0.00	1.00	negative regulation of epithelial cell proliferation involved in lung morphogenesis	BP
21	GO:0019953	2.00	Down	0.01	1.00	sexual reproduction	BP
22	GO:0006640	3.00	Up	0.01	1.00	monoacylglycerol biosynthetic process	BP
23	GO:0007507	20.00	Down	0.01	1.00	heart development	BP
24	GO:0097009	8.00	Down	0.01	1.00	energy homeostasis	BP
25	GO:0002882	1.00	Down	0.01	1.00	positive regulation of chronic inflammatory response to non-antigenic stimulus	BP
26	GO:0010753	1.00	Down	0.01	1.00	positive regulation of cGMP-mediated signaling	BP
27	GO:0043306	1.00	Down	0.01	1.00	positive regulation of mast cell degranulation	BP
28	GO:0045347	1.00	Down	0.01	1.00	negative regulation of MHC class II biosynthetic process	BP
29	GO:0044281	1.00	Down	0.01	1.00	small molecule metabolic process	BP
30	GO:0016035	1.00	Up	0.01	1.00	zeta DNA polymerase complex	CC
31	GO:0042772	1.00	Up	0.01	1.00	DNA damage response, signal transduction resulting in transcription	BP
32	GO:0070828	2.00	Down	0.01	1.00	heterochromatin organization	BP
33	GO:0034383	5.00	Up	0.01	1.00	low-density lipoprotein particle clearance	BP
34	GO:0031214	9.00	Down	0.01	1.00	biomineral tissue development	BP
35	GO:0019934	4.00	Down	0.01	1.00	cGMP-mediated signaling	BP
36	GO:0004622	6.00	Up	0.01	1.00	lysophospholipase activity	MF
37	GO:0005879	15.00	Up	0.01	1.00	axonemal microtubule	CC
38	GO:0006071	2.00	Up	0.01	1.00	glycerol metabolic process	BP
39	GO:0071404	7.00	Up	0.01	1.00	cellular response to low-density lipoprotein particle stimulus	BP
40	GO:0004649	2.00	Down	0.02	1.00	poly(ADP-ribose) glycohydrolase activity	MF
41	GO:0043405	3.00	Down	0.02	1.00	regulation of MAP kinase activity	BP
42	GO:0005501	1.00	Up	0.02	1.00	retinoid binding	MF
43	GO:0070084	1.00	Down	0.02	1.00	protein initiator methionine removal	BP
44	GO:0010572	1.00	Up	0.02	1.00	positive regulation of platelet activation	BP
45	GO:0034478	1.00	Up	0.02	1.00	phosphatidylglycerol catabolic process	BP
46	GO:0006816	10.00	Down	0.02	1.00	calcium ion transport	BP
47	GO:0008409	1.00	Up	0.02	1.00	5’-3’ exonuclease activity	MF
48	GO:0034353	1.00	Up	0.02	1.00	mRNA 5’-diphosphatase activity	MF
49	GO:0050779	1.00	Up	0.02	1.00	RNA destabilization	BP
50	GO:0071028	1.00	Up	0.02	1.00	nuclear mRNA surveillance	BP
51	GO:0097345	5.00	Down	0.02	1.00	mitochondrial outer membrane permeabilization	BP
52	GO:0051260	27.00	Up	0.02	1.00	protein homooligomerization	BP
53	GO:2000253	1.00	Up	0.02	1.00	positive regulation of feeding behavior	BP
54	GO:0019834	1.00	Down	0.02	1.00	phospholipase A2 inhibitor activity	MF
55	GO:0034115	3.00	Up	0.02	1.00	negative regulation of heterotypic cell-cell adhesion	BP
56	GO:0072557	1.00	Up	0.02	1.00	IPAF inflammasome complex	CC
57	GO:0004693	4.00	Down	0.02	1.00	cyclin-dependent protein serine/threonine kinase activity	MF
58	GO:0001833	3.00	Down	0.02	1.00	inner cell mass cell proliferation	BP
59	GO:0048863	7.00	Up	0.02	1.00	stem cell differentiation	BP
60	GO:0051651	1.00	Down	0.02	1.00	maintenance of location in cell	BP
61	GO:1905379	1.00	Down	0.02	1.00	beta-N-acetylhexosaminidase complex	CC
62	GO:0045836	1.00	Up	0.02	1.00	positive regulation of meiotic nuclear division	BP
63	GO:0070201	1.00	Up	0.02	1.00	regulation of establishment of protein localization	BP
64	GO:0042420	1.00	Down	0.02	1.00	dopamine catabolic process	BP
65	GO:0070873	2.00	Up	0.02	1.00	regulation of glycogen metabolic process	BP
66	GO:0090267	2.00	Up	0.02	1.00	positive regulation of mitotic cell cycle spindle assembly checkpoint	BP
67	GO:0006820	1.00	Down	0.02	1.00	anion transport	BP
68	GO:0003846	1.00	Up	0.02	1.00	2-acylglycerol O-acyltransferase activity	MF
69	GO:0035336	1.00	Up	0.02	1.00	long-chain fatty-acyl-CoA metabolic process	BP
70	GO:0071400	1.00	Up	0.02	1.00	cellular response to oleic acid	BP
71	GO:0097006	1.00	Up	0.02	1.00	regulation of plasma lipoprotein particle levels	BP
72	GO:0034638	4.00	Up	0.02	1.00	phosphatidylcholine catabolic process	BP
73	GO:0046983	22.00	Down	0.02	1.00	protein dimerization activity	MF
74	GO:0046339	3.00	Up	0.03	1.00	diacylglycerol metabolic process	BP
75	GO:0030521	2.00	Down	0.03	1.00	androgen receptor signaling pathway	BP
76	GO:0008308	2.00	Down	0.03	1.00	voltage-gated anion channel activity	MF
77	GO:1903146	2.00	Down	0.03	1.00	regulation of autophagy of mitochondrion	BP
78	GO:0004930	14.00	Up	0.03	1.00	G protein-coupled receptor activity	MF
79	GO:0055089	4.00	Up	0.03	1.00	fatty acid homeostasis	BP
80	GO:0004144	2.00	Up	0.03	1.00	diacylglycerol O-acyltransferase activity	MF
81	GO:0019432	2.00	Up	0.03	1.00	triglyceride biosynthetic process	BP
82	GO:0050252	2.00	Up	0.03	1.00	retinol O-fatty-acyltransferase activity	MF
83	GO:0050435	3.00	Up	0.03	1.00	amyloid-beta metabolic process	BP
84	GO:0106074	2.00	Down	0.03	1.00	aminoacyl-tRNA metabolism involved in translational fidelity	BP
85	GO:0032693	1.00	Down	0.03	1.00	negative regulation of interleukin-10 production	BP
86	GO:0017060	1.00	Up	0.03	1.00	3-galactosyl-N-acetylglucosaminide 4-alpha-L-fucosyltransferase activity	MF
87	GO:0050681	6.00	Down	0.03	1.00	nuclear androgen receptor binding	MF
88	GO:0043849	1.00	Down	0.03	1.00	Ras palmitoyltransferase activity	MF
89	GO:0004999	1.00	Down	0.03	1.00	vasoactive intestinal polypeptide receptor activity	MF
90	GO:0007260	2.00	Down	0.03	1.00	tyrosine phosphorylation of STAT protein	BP
91	GO:0031430	1.00	Down	0.03	1.00	M band	CC
92	GO:0030878	5.00	Down	0.03	1.00	thyroid gland development	BP
93	GO:0048664	2.00	Down	0.03	1.00	neuron fate determination	BP
94	GO:0004667	2.00	Up	0.03	1.00	prostaglandin-D synthase activity	MF
95	GO:0140007	1.00	Up	0.03	1.00	KICSTOR complex	CC
96	GO:0001756	6.00	Up	0.03	1.00	somitogenesis	BP
97	GO:0001837	6.00	Down	0.03	1.00	epithelial to mesenchymal transition	BP
98	GO:0042752	13.00	Down	0.03	1.00	regulation of circadian rhythm	BP
99	GO:0051028	20.00	Down	0.03	1.00	mRNA transport	BP
100	GO:0017083	2.00	Up	0.03	1.00	4-galactosyl-N-acetylglucosaminide 3-alpha-L-fucosyltransferase activity	MF
101	GO:0036505	1.00	Up	0.03	1.00	prosaposin receptor activity	MF
102	GO:0000964	1.00	Up	0.04	1.00	mitochondrial RNA 5’-end processing	BP
103	GO:0090646	1.00	Up	0.04	1.00	mitochondrial tRNA processing	BP
104	GO:2001256	3.00	Up	0.04	1.00	regulation of store-operated calcium entry	BP
105	GO:0000281	12.00	Up	0.04	1.00	mitotic cytokinesis	BP
106	GO:0004017	3.00	Up	0.04	1.00	adenylate kinase activity	MF
107	GO:0071455	2.00	Down	0.04	1.00	cellular response to hyperoxia	BP
108	GO:0006691	1.00	Up	0.04	1.00	leukotriene metabolic process	BP
109	GO:0004553	1.00	Down	0.04	1.00	hydrolase activity, hydrolyzing O-glycosyl compounds	MF
110	GO:0140290	1.00	Down	0.04	1.00	peptidyl-serine ADP-deribosylation	BP
111	GO:0140292	1.00	Down	0.04	1.00	ADP-ribosylserine hydrolase activity	MF
112	GO:0043137	1.00	Up	0.04	1.00	DNA replication, removal of RNA primer	BP
113	GO:0045722	2.00	Down	0.04	1.00	positive regulation of gluconeogenesis	BP
114	GO:1902993	1.00	Down	0.04	1.00	positive regulation of amyloid precursor protein catabolic process	BP
115	GO:0000244	3.00	Down	0.04	1.00	spliceosomal tri-snRNP complex assembly	BP
116	GO:0010762	1.00	Up	0.04	1.00	regulation of fibroblast migration	BP
117	GO:0031932	1.00	Up	0.04	1.00	TORC2 complex	CC
118	GO:0038203	1.00	Up	0.04	1.00	TORC2 signaling	BP
119	GO:0042382	3.00	Up	0.04	1.00	paraspeckles	CC
120	GO:0043715	1.00	Down	0.04	1.00	2,3-diketo-5-methylthiopentyl-1-phosphate enolase activity	MF
121	GO:0043716	1.00	Down	0.04	1.00	2-hydroxy-3-keto-5-methylthiopentenyl-1-phosphate phosphatase activity	MF
122	GO:0043874	1.00	Down	0.04	1.00	acireductone synthase activity	MF
123	GO:0034774	1.00	Up	0.04	1.00	secretory granule lumen	CC
124	GO:0090277	1.00	Up	0.04	1.00	positive regulation of peptide hormone secretion	BP
125	GO:1900738	1.00	Up	0.04	1.00	positive regulation of phospholipase C-activating G protein-coupled receptor signaling pathway	BP
126	GO:1903817	1.00	Up	0.04	1.00	negative regulation of voltage-gated potassium channel activity	BP
127	GO:1905151	1.00	Up	0.04	1.00	negative regulation of voltage-gated sodium channel activity	BP
128	GO:0034394	8.00	Down	0.04	1.00	protein localization to cell surface	BP
129	GO:0004534	1.00	Down	0.04	1.00	5’-3’ exoribonuclease activity	MF
130	GO:0031124	1.00	Down	0.04	1.00	mRNA 3’-end processing	BP
131	GO:0070932	2.00	Up	0.04	1.00	histone H3 deacetylation	BP
132	GO:0106008	1.00	Down	0.04	1.00	2-oxoglutaramate amidase activity	MF
133	GO:1990459	3.00	Up	0.04	1.00	transferrin receptor binding	MF
134	GO:0031685	1.00	Down	0.04	1.00	adenosine receptor binding	MF
135	GO:0030155	9.00	Down	0.04	1.00	regulation of cell adhesion	BP
136	GO:0070374	27.00	Down	0.04	1.00	positive regulation of ERK1 and ERK2 cascade	BP
137	GO:0010745	1.00	Up	0.04	1.00	negative regulation of macrophage derived foam cell differentiation	BP
138	GO:0030853	1.00	Up	0.04	1.00	negative regulation of granulocyte differentiation	BP
139	GO:0033034	1.00	Up	0.04	1.00	positive regulation of myeloid cell apoptotic process	BP
140	GO:0033691	1.00	Up	0.04	1.00	sialic acid binding	MF
141	GO:0045650	1.00	Up	0.04	1.00	negative regulation of macrophage differentiation	BP
142	GO:0045923	1.00	Up	0.04	1.00	positive regulation of fatty acid metabolic process	BP
143	GO:0050805	1.00	Up	0.04	1.00	negative regulation of synaptic transmission	BP
144	GO:0090317	1.00	Up	0.04	1.00	negative regulation of intracellular protein transport	BP
145	GO:0110113	1.00	Up	0.04	1.00	positive regulation of lipid transporter activity	BP
146	GO:2000279	1.00	Up	0.04	1.00	negative regulation of DNA biosynthetic process	BP
147	GO:2000467	1.00	Up	0.04	1.00	positive regulation of glycogen (starch) synthase activity	BP
148	GO:2000478	1.00	Up	0.04	1.00	positive regulation of metanephric podocyte development	BP
149	GO:2000534	1.00	Up	0.04	1.00	positive regulation of renal albumin absorption	BP
150	GO:2000584	1.00	Up	0.04	1.00	negative regulation of platelet-derived growth factor receptor-alpha signaling pathway	BP
151	GO:2000590	1.00	Up	0.04	1.00	negative regulation of metanephric mesenchymal cell migration	BP
152	GO:0004046	1.00	Down	0.04	1.00	aminoacylase activity	MF
153	GO:0016811	1.00	Down	0.04	1.00	hydrolase activity, acting on carbon-nitrogen (but not peptide) bonds, in linear amides	MF
154	GO:0043604	1.00	Down	0.04	1.00	amide biosynthetic process	BP
155	GO:0043605	1.00	Down	0.04	1.00	cellular amide catabolic process	BP
156	GO:2000275	1.00	Down	0.04	1.00	regulation of oxidative phosphorylation uncoupler activity	BP
157	GO:0010826	2.00	Up	0.04	1.00	negative regulation of centrosome duplication	BP
158	GO:0005666	1.00	Up	0.04	1.00	RNA polymerase III complex	CC
159	GO:0034362	6.00	Up	0.04	1.00	low-density lipoprotein particle	CC
160	GO:0051131	7.00	Down	0.04	1.00	chaperone-mediated protein complex assembly	BP
161	GO:0000123	4.00	Down	0.04	1.00	histone acetyltransferase complex	CC
162	GO:0065010	2.00	Down	0.04	1.00	extracellular membrane-bounded organelle	CC
163	GO:0021615	1.00	Down	0.05	1.00	glossopharyngeal nerve morphogenesis	BP
164	GO:0048645	1.00	Down	0.05	1.00	animal organ formation	BP
165	GO:0071837	1.00	Down	0.05	1.00	HMG box domain binding	MF
166	GO:1903432	1.00	Down	0.05	1.00	regulation of TORC1 signaling	BP
167	GO:0044793	2.00	Up	0.05	1.00	negative regulation by host of viral process	BP
168	GO:0036211	9.00	Down	0.05	1.00	protein modification process	BP
169	GO:1905457	1.00	Up	0.05	1.00	negative regulation of lymphoid progenitor cell differentiation	BP
170	GO:0043981	5.00	Down	0.05	1.00	histone H4-K5 acetylation	BP
171	GO:0043982	5.00	Down	0.05	1.00	histone H4-K8 acetylation	BP
172	GO:0006024	2.00	Down	0.05	1.00	glycosaminoglycan biosynthetic process	BP
173	GO:0036477	2.00	Up	0.05	1.00	somatodendritic compartment	CC
174	GO:1903942	2.00	Up	0.05	1.00	positive regulation of respiratory gaseous exchange	BP
175	GO:0050830	11.00	Up	0.05	1.00	defense response to Gram-positive bacterium	BP
176	GO:0010524	2.00	Down	0.05	1.00	positive regulation of calcium ion transport into cytosol	BP
177	GO:0060219	1.00	Up	0.05	1.00	camera-type eye photoreceptor cell differentiation	BP
178	GO:0061298	1.00	Up	0.05	1.00	retina vasculature development in camera-type eye	BP
179	GO:0045244	1.00	Down	0.05	1.00	succinate-CoA ligase complex (GDP-forming)	CC
180	GO:0034141	3.00	Down	0.05	1.00	positive regulation of toll-like receptor 3 signaling pathway	BP
181	GO:0061564	2.00	Down	0.05	1.00	axon development	BP

Session Info

sessionInfo()

EdgeR w/o Outlier

Load data

library(edgeR)
library(statmod)
library(tximport)
library(data.table)
library(tidyverse)
library(ggrepel)
library(ascii)
options(asciiType="org", width=80)
meta.data <- data.frame(library=c("SLICK_2041", "SLICK_2128", "SLICK_2129",
                                  "SLICK_2171", "SLICK_2202", "SLICK_2287",
                                  "SLICK_466", "WT_2058", "WT_2068", "WT_2089",
                                  "WT_2116", "WT_2226", "WT_2273"),
                        group=factor(rep(c("SL", "WT"), each=7, length.out = 13),
                                     levels=c("WT", "SL")))

files <- file.path('2-salmon', meta.data$library,'quant.sf')
names(files) <- meta.data$library


design <- model.matrix(~ group, data=meta.data)
colnames(design) <- make.names(gsub("group", "", colnames(design)))

## Create tx2gene from gff
txdb <- fread("0-ref/ncbi_dataset/data/GCF_002263795.2/genomic.gff",
              skip=8, fill=T, header=F, sep="\t")[V3=='mRNA',] %>%
  dplyr::select('Attr'=V9) %>%
  as.data.frame %>%
  tibble::rowid_to_column(var='rowid') %>%
  separate_rows(Attr, sep=';') %>%
  separate(Attr, into=c('key', 'value'), sep='=') %>%
  spread(key, value) %>%
  mutate(GeneID = sub(".*GeneID:([^,]*).*", "\\1", Dbxref))

tx2gene <- dplyr::select(txdb, Name, GeneID)

txi <- tximport(files, type = "salmon", tx2gene = tx2gene)

Normalize counts, Filter Genes, and Estimate Dispersions

cts <- txi$counts
normMat <- txi$length

## Obtaining per-observation scaling factors for length, adjusted to avoid
## changing the magnitude of the counts.
normMat <- normMat/exp(rowMeans(log(normMat)))
normCts <- cts/normMat

 # Computing effective library sizes from scaled counts, to account for
# composition biases between samples.
eff.lib <- calcNormFactors(normCts) * colSums(normCts)

# Combining effective library sizes with the length factors, and calculating
# offsets for a log-link GLM.
normMat <- sweep(normMat, 2, eff.lib, "*")
normMat <- log(normMat)

# Creating a DGEList object for use in edgeR.
d <- DGEList(cts) %>% scaleOffset(normMat)

# filtering using the design information
keep <- filterByExpr(d, design)

y <- d[keep, ,keep.lib.sizes=F ] %>%
  calcNormFactors %>%
  estimateDisp(design, robust=T)

unique genes

as.data.frame(txi$counts) %>%
  rownames_to_column('GeneID') %>%
gather(-GeneID, key="library", value="count") %>%
  left_join(meta.data, by='library') %>%
  filter(count > 0) %>%
  group_by(GeneID, group) %>%
  count %>%
  spread(group, n) %>%
  filter(is.na(WT) | is.na(SL)) %>%
  gather(-GeneID, key="Type", value="Samples") %>%
  filter(Samples >= 4) %>%
  ascii(digits=0)

	GeneID	Type	Samples
1	112441545	WT	6
2	112447300	WT	5
3	505979	WT	4
4	505993	WT	4
5	523830	WT	4
6	535831	WT	4
7	615310	WT	4
8	786611	WT	4
9	101902675	SL	4
10	101905578	SL	4
11	101906363	SL	4
12	104968609	SL	5
13	107132566	SL	6
14	112449108	SL	4
15	512850	SL	4
16	516736	SL	5
17	524709	SL	4
18	526745	SL	4
19	540297	SL	4
20	613925	SL	5
21	616105	SL	4
22	780786	SL	4
23	783189	SL	4
24	790886	SL	4

Unique Gene counts

as.data.frame(txi$counts) %>%
  rownames_to_column('gene') %>%
  gather(-gene, key="library", value="count") %>%
  left_join(meta.data) %>%
  filter(count > 0) %>%
  group_by(gene, group) %>%
  count %>%
  spread(group, n) %>%
  filter(is.na(WT) | is.na(SL)) %>%
  gather(-gene, key="type", value="samples") %>%
  filter(!is.na(samples)) %>%
  group_by(type, samples) %>%
  count %>%
  group_by(type) %>%
  mutate(total=sum(n)) %>%
  spread(samples, n)

Group	Total	1	2	3	4	5	6
SL	592	424	119	33	12	3	1
WT	401	307	67	19	6	1	1

Unique Gene Venn Diagram

unique.gene.venn.plot <- as.data.frame(txi$counts) %>%
  rownames_to_column('gene') %>%
gather(-gene, key="library", value="count") %>%
  left_join(meta.data) %>%
  filter(count > 0) %>%
  group_by(gene, group) %>%
  count %>%
  spread(group, n) %>%
  filter(is.na(WT) | is.na(SL)) %>%
  gather(-gene, key="type", value="count") %>%
  filter(!is.na(count)) %>%
  arrange(count) %>%
  mutate(type=fct_rev(type)) %>%
  ggplot(aes(type,1,fill=count)) +
  geom_col() +
  scale_x_discrete(breaks=c("WT", "SL"), labels=c("Wild-Type", "SLICK")) +
  scale_y_continuous(expand=c(0,0)) +
  scale_fill_distiller(palette = "PuRd") +
  guides(fill=guide_colorbar(title="Number of \nSamples")) +
  labs(y="Number of unique genes", x="Hair Coat Type") +
  theme_bw()+
  theme(
    panel.grid = element_blank(),
    axis.text = element_text(size=14),,
    axis.title = element_text(size=16),
    legend.text = element_text(size=14),
    legend.title = element_text(size=16)
    )

unique.gene.venn.plot

Filtered/Surviving gene summary
```
table(keep)
    
```

MDS

mds.data <- plotMDS(y, plot=F)

mds.plot <- cbind(meta.data, x=mds.data$x, y=mds.data$y) %>%
  mutate(library = sub('.*_', '', library)) %>%
  ggplot(aes(x,y, color=group, label=library)) +
  geom_point() +
  geom_text_repel(show.legend = F) +
  scale_color_manual(values=c( "#cc546e", "#4db598"), labels=c("Wild-Type", "SLICK")) +
  xlab(sprintf("Dimension 1 (%0.2f%%)", mds.data$var.explained[1]*100)) +
  ylab(sprintf("Dimension 2 (%0.2f%%)", mds.data$var.explained[2]*100)) +
  labs(color="Hair Coat Type") +
  theme_bw() +
  theme(legend.position = c(0.8,0.8),
        legend.background = element_rect(fill='white', color="grey"),
        panel.grid = element_blank()) +
    theme( legend.position = c(0.8,0.9),
    panel.grid.minor = element_blank(),
    axis.text = element_text(size=14),,
    axis.title = element_text(size=16),
    legend.text = element_text(size=14),
    legend.title = element_text(size=16)
    )
mds.plot

Panel MDS/Venn

library(cowplot)

plot_grid(mds.plot, unique.gene.venn.plot, labels=c("A", "B"), rel_widths = c(3,2))

Mean-Difference (MD) plots

par(mfrow = c(3,5))
for(x in 1:13){
  plotMD(y, column=x, ylim=c(-8,8))
  abline(h=0, col="red", lty=2, lwd=2)
}

Biological Coefficient of Variation (BCV)
```
plotBCV(y)
    
```

DEG LRT test

fit <- glmFit(y,design,robust=T)
lrt <- glmLRT(fit)

desc <- txdb %>%
  dplyr::select(GeneID, gene, product) %>%
  distinct(GeneID, .keep_all=T)

topTags(lrt, n=Inf) %>%
  as.data.frame %>%
  filter(PValue <= 0.01 & abs(logFC) > 0.6) %>%
  tibble::rownames_to_column("GeneID") %>%
  left_join(desc, by="GeneID")  %>%
  ascii(format=c('s', 'f', 'f','f', 'g', 'g', 's', 's'), digits=4)

	GeneID	logFC	logCPM	LR	PValue	FDR	gene	product
1	513329	2.2612	1.3791	31.0589	2.503e-08	0.0003398	LOC513329	major allergen Equ c 1%2C transcript variant X2
2	100139092	1.2267	0.5330	15.1984	9.679e-05	0.5212	GPR34	G protein-coupled receptor 34%2C transcript variant X2
3	618181	1.9566	-0.9814	14.8696	0.0001152	0.5212	GRIN2D	glutamate ionotropic receptor NMDA type subunit 2D
4	281867	-1.7706	-0.8097	14.0097	0.0001819	0.6152	INHBA	inhibin subunit beta A%2C transcript variant X2
5	511941	0.8150	3.2598	13.5964	0.0002266	0.6152	CEACAM19	carcinoembryonic antigen related cell adhesion molecule 19
6	280839	1.4570	0.5286	12.6274	0.0003801	0.8599	LHB	luteinizing hormone beta polypeptide%2C transcript variant X1
7	280863	-1.1996	1.5895	12.2083	0.0004758	0.9225	MOG	myelin oligodendrocyte glycoprotein%2C transcript variant X4
8	519145	1.4944	-0.4487	11.5354	0.0006828	1	LOC519145	beta-galactosidase-1-like protein 2%2C transcript variant X1
9	784038	-5.6681	5.2023	11.3173	0.0007679	1	LOC784038	testis-specific H1 histone
10	107131164	0.6224	3.9435	11.2429	0.0007993	1	ADAMTSL5	ADAMTS like 5%2C transcript variant X4
11	538700	-2.2637	-0.0336	10.5714	0.001149	1	BOLA-DRB2	histocompatibility complex%2C class II%2C DR beta 2%2C transcript variant X2
12	319095	-0.9615	4.1792	10.2976	0.001332	1	ADCYAP1R1	ADCYAP receptor type I%2C transcript variant X5
13	783257	-0.6324	2.6277	10.1605	0.001435	1	ABCB6	ATP binding cassette subfamily B member 6
14	100336825	0.9281	1.0589	9.8293	0.001718	1	TTLL10	tubulin tyrosine ligase like 10%2C transcript variant X2
15	785756	-2.7780	3.8089	9.6565	0.001887	1	LOC785756	androgen binding protein beta-like
16	529491	-0.6692	2.3837	9.5262	0.002026	1	PMS1	PMS1 homolog 1%2C mismatch repair system component
17	782348	1.3334	-0.2756	9.3383	0.002244	1	LOC782348	uncharacterized LOC782348%2C transcript variant X1
18	767895	0.7918	2.5027	9.2474	0.002358	1	WDR60	WD repeat domain 60%2C transcript variant X4
19	286791	-2.0579	1.8573	9.1437	0.002496	1	BDA20	major allergen BDA20%2C transcript variant X1
20	788425	1.3382	-0.1399	8.7300	0.00313	1	LOC788425	aflatoxin B1 aldehyde reductase member 4%2C transcript variant X2
21	100298265	1.5069	0.0679	8.6451	0.003279	1	CPAMD8	C3 and PZP like%2C alpha-2-macroglobulin domain containing 8%2C transcript variant X1
22	518047	-2.0065	-0.2084	8.6408	0.003287	1	SLCO1C1	solute carrier organic anion transporter family member 1C1%2C transcript variant X2
23	517459	-3.0599	2.1857	8.5696	0.003418	1	GABRP	gamma-aminobutyric acid type A receptor pi subunit%2C transcript variant X2
24	615323	0.8781	5.2371	8.5691	0.003419	1	GRP	gastrin releasing peptide
25	539733	0.8105	2.2088	8.4830	0.003585	1	RWDD2A	RWD domain containing 2A%2C transcript variant X1
26	534505	-1.9690	1.8282	8.4545	0.003641	1	CPXM2	carboxypeptidase X%2C M14 family member 2
27	505551	1.1909	-0.8951	8.3957	0.003761	1	MYCBPAP	MYCBP associated protein
28	519940	-1.6339	1.5724	8.3492	0.003859	1	RDH16	retinol dehydrogenase 16 (all-trans)%2C transcript variant X5
29	539467	0.9590	3.2906	8.2536	0.004067	1	PKDCC	protein kinase domain containing%2C cytoplasmic
30	100271839	-0.7476	1.3874	7.6771	0.005593	1	C23H6orf141	chromosome 23 C6orf141 homolog
31	530346	1.1875	0.0604	7.6219	0.005766	1	ZNF621	zinc finger protein 621
32	286858	1.1112	8.5204	7.6061	0.005817	1	PTGDS	prostaglandin D2 synthase%2C transcript variant X2
33	511042	-1.1189	1.7152	7.5741	0.005921	1	NMNAT2	nicotinamide nucleotide adenylyltransferase 2%2C transcript variant X4
34	615896	-2.3026	1.4550	7.4777	0.006247	1	SYCP3	synaptonemal complex protein 3%2C transcript variant X2
35	505830	0.9327	0.3466	7.3412	0.006739	1	APOM	apolipoprotein M%2C transcript variant X2
36	280948	1.0145	3.9968	7.1634	0.00744	1	TTR	transthyretin
37	524150	-1.5792	-1.6256	6.8693	0.008769	1	HOXA2	homeobox A2
38	100335365	1.2305	-1.3344	6.8628	0.008801	1	STPG4	sperm-tail PG-rich repeat containing 4%2C transcript variant X2
39	404129	1.4760	-0.5501	6.8321	0.008953	1	DGAT2	diacylglycerol O-acyltransferase 2
40	100848815	-1.4578	5.5645	6.7733	0.009253	1	LOC100848815	SLA class II histocompatibility antigen%2C DQ haplotype D alpha chain-like
41	534578	0.6853	3.6240	6.7374	0.009441	1	LOC534578	vascular cell adhesion molecule 1-like
42	507280	0.9027	4.2938	6.7097	0.009589	1	TM4SF18	transmembrane 4 L six family member 18%2C transcript variant X1
43	613448	1.2015	0.9361	6.6744	0.009781	1	AK5	adenylate kinase 5

MD plot

names <- dplyr::select(txdb, GeneID, gene) %>%
  distinct(GeneID, .keep_all=T)

topTags(lrt, n=Inf)$table %>%
    as.data.frame %>%
    rownames_to_column("GeneID") %>%
    left_join(names) %>%
    mutate(Sig = (PValue <= 0.01 & abs(logFC) >= 0.6),
          Dir = ifelse(Sig, ifelse(logFC < 0, "Down", "Up"), "Not Significant"),
          label=if_else(Sig, gene, NA)) %>%
     arrange(Sig) %>%
ggplot(aes(logCPM, logFC, color=Dir, size=Sig, label=label)) +
  geom_point() +
  geom_label_repel(size = 4, max.overlaps=100) +
  scale_size_manual(values=c(2/5, 1), guide=F) +
  scale_color_manual(values=c(Up="red", Down="blue", "grey"),
                     labels=c(Up="Up (25)", Down="Down (12)")) +
  xlab(expression("log"[2]*"CPM")) +
  ylab(expression("log"[2]*"FC")) +
  labs(color="Significant\n(p-value <= 0.01, |logFC| >= 0.6)") +
  theme_bw() +
  theme( legend.position = c(0.8,0.9),
    panel.grid.minor = element_blank(),
    axis.text = element_text(size=14),,
    axis.title = element_text(size=16),
    legend.text = element_text(size=14),
    legend.title = element_text(size=16)
    )

Volcano

topTags(lrt, n=Inf)$table %>%
   mutate(Sig = (PValue <= 0.01 & abs(logFC) >= 0.6),
          Dir = ifelse(Sig, ifelse(logFC < 0, "Down", "Up"), "Not Significant")) %>%
     arrange(Sig) %>%
ggplot(aes(logFC, -log(PValue,10), color=Dir, size=Sig)) + # -log10 conversion
  geom_point() +
  scale_size_manual(values=c(2/5, 1), guide=F) +
  scale_color_manual(values=c(Up="red", Down="blue", "grey")) +
  xlab(expression("log"[2]*"FC")) +
  ylab(expression("-log"[10]*"p-value")) +
  theme_minimal() +
  theme(legend.position = 'none')

KEGG

Download pathways and run fry

kg.list <- getGeneKEGGLinks("bta")
kg.sets <- split(kg.list$GeneID, kg.list$PathwayID)

kg.names <- getKEGGPathwayNames('bta')

kegg.fry <- fry(y, kg.sets, design)

Get pathways with p-value <= 0.05

kegg.fry %>%
  rownames_to_column('PathwayID') %>%
  mutate(PathwayID=sub('path:', '', PathwayID)) %>%
  left_join(kg.names, by="PathwayID") %>%
  filter(PValue <= 0.05 | PValue.Mixed <= 0.05) %>%
  mutate(Description = sub(' - Bos taurus \\(cow\\)', '', Description)) %>%
  ascii(format=c('s', 'd', 's', 'g' ,'g' ,'g' , 'g', 's'))

	PathwayID	NGenes	Direction	PValue	FDR	PValue.Mixed	FDR.Mixed	Description
1	bta00515	21	Up	0.00057	0.2	0.31	0.93	Mannose type O-glycan biosynthesis
2	bta00562	67	Down	0.012	0.99	0.67	0.93	Inositol phosphate metabolism
3	bta04070	86	Down	0.014	0.99	0.62	0.93	Phosphatidylinositol signaling system
4	bta00590	44	Up	0.021	0.99	0.48	0.93	Arachidonic acid metabolism
5	bta04530	146	Down	0.027	0.99	0.53	0.93	Tight junction
6	bta05211	67	Down	0.041	0.99	0.72	0.93	Renal cell carcinoma
7	bta00310	57	Down	0.048	0.99	0.85	0.94	Lysine degradation
8	bta04012	72	Down	0.05	0.99	0.69	0.93	ErbB signaling pathway
9	bta05033	9	Up	0.24	0.99	0.0045	0.93	Nicotine addiction

bta00515

Mannose type O-glycan biosynthesis

library(pathview)
pathview(pathway.id=pathway,
         kegg.dir = '4-kegg/',
         gene.data = setNames(lrt$table$logFC, row.names(lrt$table)),
         species='bta')

bta00562

Inositol phosphate metabolism

bta04070

Phosphatidylinositol signaling system

bta04530

Tight junction

bta04012

ErbB signaling pathway

-

topTags(lrt, n=Inf) %>%
  as.data.frame %>%
  tibble::rownames_to_column("GeneID") %>%
  filter(GeneID %in% kg.sets[['path:bta00590']]) %>%
  left_join(desc, by="GeneID")  %>%
  ascii(format=c('s', 'f', 'f','f', 'g', 'g', 's', 's'), digits=4)

GO

library(org.Bt.eg.db)
library(GO.db)

go.db <- select(org.Bt.eg.db,
       keys=keys(org.Bt.eg.db, keytype = "ENTREZID"),
       columns="GO",
       keytype="ENTREZID") %>%
  filter(!is.na(GO))


go.sets <- split(go.db$ENTREZID, go.db$GO)

go.names <- select(GO.db,
                   keys=names(go.sets),
                   columns=c("GOID", "TERM", "ONTOLOGY"),
                   keytype="GOID")

go.fry <- fry(y, go.sets, design)

Get go terms with p-value <= 0.05

go.fry %>%
  rownames_to_column('GOID') %>%
  filter(PValue <= 0.05) %>%
  left_join(go.names, by='GOID') %>%
  dplyr::select(-ends_with('.Mixed')) %>%
  ascii(format=c('s', 'd', 's', 'g', 'g', 's', 's'))

	GOID	NGenes	Direction	PValue	FDR	TERM	ONTOLOGY
1	GO:0045648	6	Down	0.00056	1	positive regulation of erythrocyte differentiation	BP
2	GO:1903598	3	Down	0.00073	1	positive regulation of gap junction assembly	BP
3	GO:0007507	20	Down	0.0013	1	heart development	BP
4	GO:0009409	6	Down	0.0017	1	response to cold	BP
5	GO:0021773	1	Down	0.0019	1	striatal medium spiny neuron differentiation	BP
6	GO:0035987	1	Down	0.0019	1	endodermal cell differentiation	BP
7	GO:0042541	1	Down	0.0019	1	hemoglobin biosynthetic process	BP
8	GO:0042701	1	Down	0.0019	1	progesterone secretion	BP
9	GO:0043509	1	Down	0.0019	1	activin A complex	CC
10	GO:0043512	1	Down	0.0019	1	inhibin A complex	CC
11	GO:0046880	1	Down	0.0019	1	regulation of follicle-stimulating hormone secretion	BP
12	GO:0048333	1	Down	0.0019	1	mesodermal cell differentiation	BP
13	GO:0051799	1	Down	0.0019	1	negative regulation of hair follicle development	BP
14	GO:0070699	1	Down	0.0019	1	type II activin receptor binding	MF
15	GO:0097154	1	Down	0.0019	1	GABAergic neuron differentiation	BP
16	GO:0140768	2	Down	0.0024	1	protein ADP-ribosyltransferase-substrate adaptor activity	MF
17	GO:0046628	1	Down	0.0025	1	positive regulation of insulin receptor signaling pathway	BP
18	GO:0008308	2	Down	0.0027	1	voltage-gated anion channel activity	MF
19	GO:0001942	5	Down	0.0039	1	hair follicle development	BP
20	GO:0033962	3	Down	0.0045	1	P-body assembly	BP
21	GO:0005643	5	Down	0.0047	1	nuclear pore	CC
22	GO:0002882	1	Down	0.0055	1	positive regulation of chronic inflammatory response to non-antigenic stimulus	BP
23	GO:0010753	1	Down	0.0055	1	positive regulation of cGMP-mediated signaling	BP
24	GO:0043306	1	Down	0.0055	1	positive regulation of mast cell degranulation	BP
25	GO:0050693	2	Down	0.0061	1	LBD domain binding	MF
26	GO:0051302	6	Up	0.0066	1	regulation of cell division	BP
27	GO:0007616	5	Up	0.0069	1	long-term memory	BP
28	GO:0034968	2	Up	0.0073	1	histone lysine methylation	BP
29	GO:0030151	2	Down	0.0075	1	molybdenum ion binding	MF
30	GO:0032793	2	Up	0.0079	1	positive regulation of CREB transcription factor activity	BP
31	GO:0008265	1	Down	0.0085	1	Mo-molybdopterin cofactor sulfurase activity	MF
32	GO:0043545	1	Down	0.0085	1	molybdopterin cofactor metabolic process	BP
33	GO:0102867	1	Down	0.0085	1	molybdenum cofactor sulfurtransferase activity	MF
34	GO:0000028	1	Up	0.0085	1	ribosomal small subunit assembly	BP
35	GO:0005388	2	Down	0.0087	1	P-type calcium transporter activity	MF
36	GO:0021740	1	Down	0.0089	1	principal sensory nucleus of trigeminal nerve development	BP
37	GO:0021960	1	Down	0.0089	1	anterior commissure morphogenesis	BP
38	GO:0044300	1	Down	0.0089	1	cerebellar mossy fiber	CC
39	GO:0060486	1	Down	0.0089	1	club cell differentiation	BP
40	GO:0060509	1	Down	0.0089	1	type I pneumocyte differentiation	BP
41	GO:0061141	1	Down	0.0089	1	lung ciliated cell differentiation	BP
42	GO:0071679	1	Down	0.0089	1	commissural neuron axon guidance	BP
43	GO:2000791	1	Down	0.0089	1	negative regulation of mesenchymal cell proliferation involved in lung development	BP
44	GO:2000795	1	Down	0.0089	1	negative regulation of epithelial cell proliferation involved in lung morphogenesis	BP
45	GO:0019953	2	Down	0.009	1	sexual reproduction	BP
46	GO:0048664	2	Down	0.0093	1	neuron fate determination	BP
47	GO:0008242	1	Up	0.0096	1	omega peptidase activity	MF
48	GO:0070084	1	Down	0.011	1	protein initiator methionine removal	BP
49	GO:0006640	3	Up	0.011	1	monoacylglycerol biosynthetic process	BP
50	GO:0001754	1	Down	0.011	1	eye photoreceptor cell differentiation	BP
51	GO:0033210	1	Down	0.011	1	leptin-mediated signaling pathway	BP
52	GO:0042755	1	Down	0.011	1	eating behavior	BP
53	GO:0044320	1	Down	0.011	1	cellular response to leptin stimulus	BP
54	GO:0048708	1	Down	0.011	1	astrocyte differentiation	BP
55	GO:0060019	1	Down	0.011	1	radial glial cell differentiation	BP
56	GO:0070102	1	Down	0.011	1	interleukin-6-mediated signaling pathway	BP
57	GO:0017060	1	Up	0.012	1	3-galactosyl-N-acetylglucosaminide 4-alpha-L-fucosyltransferase activity	MF
58	GO:0097009	8	Down	0.012	1	energy homeostasis	BP
59	GO:2000637	5	Down	0.012	1	positive regulation of miRNA-mediated gene silencing	BP
60	GO:0051260	29	Up	0.013	1	protein homooligomerization	BP
61	GO:0042420	1	Down	0.014	1	dopamine catabolic process	BP
62	GO:0017083	2	Up	0.014	1	4-galactosyl-N-acetylglucosaminide 3-alpha-L-fucosyltransferase activity	MF
63	GO:0008409	1	Up	0.014	1	5’-3’ exonuclease activity	MF
64	GO:0034353	1	Up	0.014	1	mRNA 5’-diphosphatase activity	MF
65	GO:0050779	1	Up	0.014	1	RNA destabilization	BP
66	GO:0071028	1	Up	0.014	1	nuclear mRNA surveillance	BP
67	GO:0034774	1	Up	0.014	1	secretory granule lumen	CC
68	GO:0090277	1	Up	0.014	1	positive regulation of peptide hormone secretion	BP
69	GO:1900738	1	Up	0.014	1	positive regulation of phospholipase C-activating G protein-coupled receptor signaling pathway	BP
70	GO:1903817	1	Up	0.014	1	negative regulation of voltage-gated potassium channel activity	BP
71	GO:1905151	1	Up	0.014	1	negative regulation of voltage-gated sodium channel activity	BP
72	GO:0007041	3	Down	0.015	1	lysosomal transport	BP
73	GO:0044281	1	Down	0.015	1	small molecule metabolic process	BP
74	GO:0070934	3	Down	0.016	1	CRD-mediated mRNA stabilization	BP
75	GO:0008526	2	Up	0.016	1	phosphatidylinositol transfer activity	MF
76	GO:0120019	2	Up	0.016	1	phosphatidylcholine transfer activity	MF
77	GO:0070828	2	Down	0.017	1	heterochromatin organization	BP
78	GO:0019934	4	Down	0.017	1	cGMP-mediated signaling	BP
79	GO:0006357	39	Down	0.017	1	regulation of transcription by RNA polymerase II	BP
80	GO:0045722	2	Down	0.017	1	positive regulation of gluconeogenesis	BP
81	GO:0045347	1	Down	0.018	1	negative regulation of MHC class II biosynthetic process	BP
82	GO:0004999	1	Down	0.018	1	vasoactive intestinal polypeptide receptor activity	MF
83	GO:0055089	4	Up	0.018	1	fatty acid homeostasis	BP
84	GO:0051452	1	Down	0.019	1	intracellular pH reduction	BP
85	GO:0000086	6	Up	0.019	1	G2/M transition of mitotic cell cycle	BP
86	GO:0016035	1	Up	0.02	1	zeta DNA polymerase complex	CC
87	GO:0042772	1	Up	0.02	1	DNA damage response, signal transduction resulting in transcription	BP
88	GO:0071795	1	Down	0.02	1	K11-linked polyubiquitin modification-dependent protein binding	MF
89	GO:0043405	3	Down	0.02	1	regulation of MAP kinase activity	BP
90	GO:0071455	2	Down	0.02	1	cellular response to hyperoxia	BP
91	GO:1903942	2	Up	0.02	1	positive regulation of respiratory gaseous exchange	BP
92	GO:0034383	5	Up	0.02	1	low-density lipoprotein particle clearance	BP
93	GO:0031214	9	Down	0.022	1	biomineral tissue development	BP
94	GO:0000281	12	Up	0.022	1	mitotic cytokinesis	BP
95	GO:0060219	1	Up	0.023	1	camera-type eye photoreceptor cell differentiation	BP
96	GO:0061298	1	Up	0.023	1	retina vasculature development in camera-type eye	BP
97	GO:0031781	1	Up	0.023	1	type 3 melanocortin receptor binding	MF
98	GO:0031782	1	Up	0.023	1	type 4 melanocortin receptor binding	MF
99	GO:0032098	1	Up	0.023	1	regulation of appetite	BP
100	GO:0033059	1	Up	0.023	1	cellular pigmentation	BP
101	GO:0070996	1	Up	0.023	1	type 1 melanocortin receptor binding	MF
102	GO:0140668	1	Up	0.023	1	positive regulation of oxytocin production	BP
103	GO:1990680	1	Up	0.023	1	response to melanocyte-stimulating hormone	BP
104	GO:0106008	1	Down	0.023	1	2-oxoglutaramate amidase activity	MF
105	GO:1903670	2	Down	0.024	1	regulation of sprouting angiogenesis	BP
106	GO:0045064	2	Up	0.024	1	T-helper 2 cell differentiation	BP
107	GO:0043011	3	Up	0.024	1	myeloid dendritic cell differentiation	BP
108	GO:1903546	2	Up	0.025	1	protein localization to photoreceptor outer segment	BP
109	GO:0006820	1	Down	0.025	1	anion transport	BP
110	GO:0042382	3	Up	0.025	1	paraspeckles	CC
111	GO:0006071	2	Up	0.026	1	glycerol metabolic process	BP
112	GO:0004622	6	Up	0.026	1	lysophospholipase activity	MF
113	GO:0070873	2	Up	0.026	1	regulation of glycogen metabolic process	BP
114	GO:0004144	2	Up	0.027	1	diacylglycerol O-acyltransferase activity	MF
115	GO:0019432	2	Up	0.027	1	triglyceride biosynthetic process	BP
116	GO:0050252	2	Up	0.027	1	retinol O-fatty-acyltransferase activity	MF
117	GO:0031685	1	Down	0.027	1	adenosine receptor binding	MF
118	GO:1990756	11	Down	0.027	1	ubiquitin ligase-substrate adaptor activity	MF
119	GO:2000253	1	Up	0.029	1	positive regulation of feeding behavior	BP
120	GO:0060397	2	Down	0.03	1	growth hormone receptor signaling pathway via JAK-STAT	BP
121	GO:0046983	22	Down	0.03	1	protein dimerization activity	MF
122	GO:0005879	15	Up	0.031	1	axonemal microtubule	CC
123	GO:1903146	2	Down	0.031	1	regulation of autophagy of mitochondrion	BP
124	GO:0071404	7	Up	0.031	1	cellular response to low-density lipoprotein particle stimulus	BP
125	GO:0004649	2	Down	0.031	1	poly(ADP-ribose) glycohydrolase activity	MF
126	GO:0003846	1	Up	0.031	1	2-acylglycerol O-acyltransferase activity	MF
127	GO:0035336	1	Up	0.031	1	long-chain fatty-acyl-CoA metabolic process	BP
128	GO:0071400	1	Up	0.031	1	cellular response to oleic acid	BP
129	GO:0097006	1	Up	0.031	1	regulation of plasma lipoprotein particle levels	BP
130	GO:0007260	2	Down	0.031	1	tyrosine phosphorylation of STAT protein	BP
131	GO:0097345	5	Down	0.032	1	mitochondrial outer membrane permeabilization	BP
132	GO:1990459	3	Up	0.032	1	transferrin receptor binding	MF
133	GO:0036143	1	Down	0.033	1	kringle domain binding	MF
134	GO:0005501	1	Up	0.033	1	retinoid binding	MF
135	GO:0004176	2	Down	0.033	1	ATP-dependent peptidase activity	MF
136	GO:0030878	5	Down	0.033	1	thyroid gland development	BP
137	GO:0006816	10	Down	0.033	1	calcium ion transport	BP
138	GO:0004693	4	Down	0.034	1	cyclin-dependent protein serine/threonine kinase activity	MF
139	GO:0044770	2	Up	0.034	1	cell cycle phase transition	BP
140	GO:0010572	1	Up	0.035	1	positive regulation of platelet activation	BP
141	GO:0034478	1	Up	0.035	1	phosphatidylglycerol catabolic process	BP
142	GO:0004550	6	Up	0.035	1	nucleoside diphosphate kinase activity	MF
143	GO:0043137	1	Up	0.035	1	DNA replication, removal of RNA primer	BP
144	GO:0008294	1	Up	0.036	1	calcium- and calmodulin-responsive adenylate cyclase activity	MF
145	GO:0099061	1	Up	0.036	1	integral component of postsynaptic density membrane	CC
146	GO:0150076	1	Up	0.036	1	neuroinflammatory response	BP
147	GO:1900273	1	Up	0.036	1	positive regulation of long-term synaptic potentiation	BP
148	GO:0046339	3	Up	0.036	1	diacylglycerol metabolic process	BP
149	GO:0031648	12	Up	0.036	1	protein destabilization	BP
150	GO:0072539	1	Up	0.036	1	T-helper 17 cell differentiation	BP
151	GO:0015936	2	Down	0.036	1	coenzyme A metabolic process	BP
152	GO:1902993	1	Down	0.036	1	positive regulation of amyloid precursor protein catabolic process	BP
153	GO:0034394	8	Down	0.037	1	protein localization to cell surface	BP
154	GO:2001241	3	Down	0.037	1	positive regulation of extrinsic apoptotic signaling pathway in absence of ligand	BP
155	GO:0048863	7	Up	0.037	1	stem cell differentiation	BP
156	GO:0090190	2	Down	0.038	1	positive regulation of branching involved in ureteric bud morphogenesis	BP
157	GO:0050908	3	Up	0.038	1	detection of light stimulus involved in visual perception	BP
158	GO:2000147	3	Up	0.038	1	positive regulation of cell motility	BP
159	GO:0006493	7	Up	0.038	1	protein O-linked glycosylation	BP
160	GO:0071425	4	Down	0.039	1	hematopoietic stem cell proliferation	BP
161	GO:0050435	3	Up	0.039	1	amyloid-beta metabolic process	BP
162	GO:0050861	2	Up	0.039	1	positive regulation of B cell receptor signaling pathway	BP
163	GO:0019834	1	Down	0.04	1	phospholipase A2 inhibitor activity	MF
164	GO:0015760	2	Up	0.04	1	glucose-6-phosphate transport	BP
165	GO:1904888	3	Down	0.04	1	cranial skeletal system development	BP
166	GO:0072557	1	Up	0.041	1	IPAF inflammasome complex	CC
167	GO:0030521	2	Down	0.041	1	androgen receptor signaling pathway	BP
168	GO:0031430	1	Down	0.041	1	M band	CC
169	GO:0120187	2	Up	0.042	1	positive regulation of protein localization to chromatin	BP
170	GO:0045836	1	Up	0.042	1	positive regulation of meiotic nuclear division	BP
171	GO:0070201	1	Up	0.042	1	regulation of establishment of protein localization	BP
172	GO:0034115	3	Up	0.042	1	negative regulation of heterotypic cell-cell adhesion	BP
173	GO:1904322	3	Up	0.042	1	cellular response to forskolin	BP
174	GO:1990357	2	Down	0.043	1	terminal web	CC
175	GO:0000118	2	Up	0.044	1	histone deacetylase complex	CC
176	GO:0090267	2	Up	0.044	1	positive regulation of mitotic cell cycle spindle assembly checkpoint	BP
177	GO:0010745	1	Up	0.044	1	negative regulation of macrophage derived foam cell differentiation	BP
178	GO:0030853	1	Up	0.044	1	negative regulation of granulocyte differentiation	BP
179	GO:0033034	1	Up	0.044	1	positive regulation of myeloid cell apoptotic process	BP
180	GO:0033691	1	Up	0.044	1	sialic acid binding	MF
181	GO:0045650	1	Up	0.044	1	negative regulation of macrophage differentiation	BP
182	GO:0045923	1	Up	0.044	1	positive regulation of fatty acid metabolic process	BP
183	GO:0050805	1	Up	0.044	1	negative regulation of synaptic transmission	BP
184	GO:0090317	1	Up	0.044	1	negative regulation of intracellular protein transport	BP
185	GO:0110113	1	Up	0.044	1	positive regulation of lipid transporter activity	BP
186	GO:2000279	1	Up	0.044	1	negative regulation of DNA biosynthetic process	BP
187	GO:2000467	1	Up	0.044	1	positive regulation of glycogen (starch) synthase activity	BP
188	GO:2000478	1	Up	0.044	1	positive regulation of metanephric podocyte development	BP
189	GO:2000534	1	Up	0.044	1	positive regulation of renal albumin absorption	BP
190	GO:2000584	1	Up	0.044	1	negative regulation of platelet-derived growth factor receptor-alpha signaling pathway	BP
191	GO:2000590	1	Up	0.044	1	negative regulation of metanephric mesenchymal cell migration	BP
192	GO:0051258	4	Up	0.045	1	protein polymerization	BP
193	GO:0001833	3	Down	0.045	1	inner cell mass cell proliferation	BP
194	GO:0051651	1	Down	0.045	1	maintenance of location in cell	BP
195	GO:1905379	1	Down	0.045	1	beta-N-acetylhexosaminidase complex	CC
196	GO:0001933	13	Up	0.046	1	negative regulation of protein phosphorylation	BP
197	GO:0060279	2	Down	0.046	1	positive regulation of ovulation	BP
198	GO:0001837	6	Down	0.047	1	epithelial to mesenchymal transition	BP
199	GO:2001256	3	Up	0.047	1	regulation of store-operated calcium entry	BP
200	GO:0043970	1	Up	0.047	1	histone H3-K9 acetylation	BP
201	GO:0050830	11	Up	0.048	1	defense response to Gram-positive bacterium	BP
202	GO:0051603	3	Up	0.048	1	proteolysis involved in protein catabolic process	BP
203	GO:0003085	2	Up	0.048	1	negative regulation of systemic arterial blood pressure	BP
204	GO:0032693	1	Down	0.048	1	negative regulation of interleukin-10 production	BP
205	GO:0034141	3	Down	0.049	1	positive regulation of toll-like receptor 3 signaling pathway	BP
206	GO:0009154	1	Up	0.049	1	purine ribonucleotide catabolic process	BP
207	GO:0016929	1	Up	0.049	1	deSUMOylase activity	MF
208	GO:0036505	1	Up	0.049	1	prosaposin receptor activity	MF
209	GO:0042310	2	Down	0.049	1	vasoconstriction	BP
210	GO:0036435	2	Down	0.049	1	K48-linked polyubiquitin modification-dependent protein binding	MF
211	GO:0065010	2	Down	0.049	1	extracellular membrane-bounded organelle	CC
212	GO:0090049	1	Down	0.05	1	regulation of cell migration involved in sprouting angiogenesis	BP
213	GO:1903026	1	Down	0.05	1	negative regulation of RNA polymerase II regulatory region sequence-specific DNA binding	BP
214	GO:1903378	1	Down	0.05	1	positive regulation of oxidative stress-induced neuron intrinsic apoptotic signaling pathway	BP
215	GO:1903955	1	Down	0.05	1	positive regulation of protein targeting to mitochondrion	BP
216	GO:1990452	1	Down	0.05	1	Parkin-FBXW7-Cul1 ubiquitin ligase complex	CC
217	GO:0043849	1	Down	0.05	1	Ras palmitoyltransferase activity	MF

GO terms of interest gene expression

GO:0007507: heart development
GO:0045648: positive regulation of erythrocyte differentiation
GO:0007507: heart development
GO:0009409: response to cold
GO:0001942: hair follicle development
GO:0097009: energy homeostasis
GO:0051260: protein homooligomerization
GO:0140668: positive regulation of oxytocin production
GO:0060397: growth hormone receptor signaling pathway via JAK-STAT
GO:0006816: calcium ion transport
GO:0050830: defense response to Gram-positive bacterium
GO:0003085: negative regulation of systemic arterial blood pressure
GO:0042310: vasoconstriction
GO:0051302: regulation of cell division

counts.per.million <- cpm(y) %>%
  as.data.frame %>%
  rownames_to_column("GeneID")

go.data <- go.terms.of.interest %>%
  separate(V1, into=c('id', 'term'), sep=" :: ") %>%
  mutate(name = paste(id, term, sep=' - ')) %>%
  pull(id, name=name) %>%
  lapply(function(x){
    filter(counts.per.million, GeneID %in% go.sets[[x]]) %>%
      left_join(desc, by="GeneID")
    })


 mapply(function(data, name) ascii(data, name=name), go.data, names(go.data), SIMPLIFY = F)


## go.plots <- mapply(function(data, name) {
##   gather(data, -GeneID, key='sample', value='cpm') %>%
##   left_join(desc, by="GeneID")  %>%
##     ggplot(aes(sample, product, fill=cpm)) +
##     geom_tile() +
##     ggtitle(name)
## }, go.data, names(go.data), SIMPLIFY = F)

	GeneID	SLICK_2041	SLICK_2128	SLICK_2129	SLICK_2171	SLICK_2202	SLICK_2287	SLICK_466	WT_2058	WT_2068	WT_2089	WT_2116	WT_2226	WT_2273	gene	product
1	280939	16.94	14.04	20.79	11.52	16.07	14.08	15.52	13.27	21.43	11.18	29.90	18.78	16.79	TEK	TEK receptor tyrosine kinase
2	281134	1.86	3.77	6.49	2.30	5.05	6.06	4.63	1.66	4.81	2.00	2.56	3.78	3.95	ECE2	endothelin converting enzyme 2%2C transcript variant 2a
3	281193	7.28	10.21	6.89	10.15	5.37	7.58	6.03	10.50	7.01	12.63	8.67	5.40	6.48	GJA1	gap junction protein alpha 1
4	281270	2.70	1.46	1.59	4.59	1.34	2.90	3.14	3.79	1.18	4.91	3.98	5.02	3.18	HOPX	HOP homeobox%2C transcript variant X3
5	281733	5.07	3.13	3.26	8.76	4.79	7.57	6.95	9.17	4.72	5.98	5.80	4.56	9.08	CXADR	CXADR%2C Ig-like cell adhesion molecule%2C transcript variant X2
6	282321	12.91	13.52	6.38	13.61	9.82	12.60	13.43	16.22	11.38	13.84	15.60	8.61	11.73	PPP3R1	protein phosphatase 3 regulatory subunit B%2C alpha
7	282382	28.25	33.40	32.41	22.80	36.04	41.62	34.63	29.67	37.23	25.04	30.00	26.82	34.54	TGFBR1	transforming growth factor beta receptor 1%2C transcript variant X2
8	282682	96.48	118.55	104.65	58.00	113.50	115.97	96.90	82.98	118.92	63.92	101.75	145.76	110.85	GRK2	G protein-coupled receptor kinase 2%2C transcript variant X1
9	338068	13.62	12.41	7.58	17.67	10.74	12.82	10.20	19.40	10.65	21.30	14.80	9.09	13.11	ACVR1	activin A receptor type 1%2C transcript variant X1
10	506085	6.25	6.10	7.41	7.26	6.38	3.96	3.94	5.50	4.90	5.65	7.26	6.91	7.01	ZDHHC16	zinc finger DHHC-type containing 16%2C transcript variant X3
11	509656	13.93	20.62	19.93	13.90	17.87	16.26	16.60	15.71	19.94	14.95	19.50	22.53	19.89	PSKH1	protein serine kinase H1
12	512254	14.28	14.09	11.33	14.39	11.81	10.63	13.29	17.31	12.14	12.17	24.27	10.09	13.61	MEF2C	myocyte enhancer factor 2C%2C transcript variant X9
13	514997	8.78	9.24	15.24	7.71	8.89	10.06	7.33	5.61	11.24	7.86	5.13	8.99	8.30	ZNHIT1	zinc finger HIT-type containing 1
14	534041	11.77	12.38	18.45	10.11	11.45	14.24	16.55	15.36	19.48	10.40	12.63	19.27	19.33	MICAL2	microtubule associated monooxygenase%2C calponin and LIM domain containing 2%2C transcript variant X1
15	534069	34.58	37.30	23.76	38.46	29.75	29.27	27.73	48.82	27.82	39.82	32.13	33.18	33.27	TGFB2	transforming growth factor beta 2%2C transcript variant X3
16	538690	11.47	14.06	21.67	11.65	19.24	5.79	17.01	9.73	20.46	12.49	15.66	15.03	15.97	ID3	inhibitor of DNA binding 3%2C HLH protein
17	539430	17.64	22.10	19.65	15.66	27.24	15.05	19.45	19.32	18.72	20.67	10.97	14.14	8.76	OSR1	odd-skipped related transciption factor 1
18	539696	21.82	21.32	20.85	21.81	26.57	23.31	24.24	20.41	28.90	22.22	27.81	21.44	22.11	HEXIM1	hexamethylene bisacetamide inducible 1
19	615492	6.96	6.87	5.74	8.43	5.27	5.14	8.09	4.82	9.53	8.39	4.00	7.62	5.34	SMG9	SMG9%2C nonsense mediated mRNA decay factor
20	786335	18.30	16.31	14.01	20.27	16.04	21.45	12.42	22.04	17.71	21.72	15.49	15.65	22.29	GYS1	glycogen synthase 1

$`GO:0045648 - positive regulation of erythrocyte differentiation`

	GeneID	SLICK_2041	SLICK_2128	SLICK_2129	SLICK_2171	SLICK_2202	SLICK_2287	SLICK_466	WT_2058	WT_2068	WT_2089	WT_2116	WT_2226	WT_2273	gene	product
1	281598	10.76	11.18	8.62	7.48	9.25	11.68	13.72	15.69	9.27	10.06	19.39	8.22	13.77	ACVR2A	activin A receptor type 2A
2	281867	0.19	0.21	0.24	0.40	0.32	0.12	0.08	0.48	1.34	1.31	0.42	0.76	0.44	INHBA	inhibin subunit beta A%2C transcript variant X2
3	281871	0.81	0.34	0.86	2.72	1.52	45.47	1.10	2.38	1.74	3.06	6.41	2.00	1.71	ISG15	ISG15 ubiquitin-like modifier
4	282376	48.07	59.84	50.23	25.17	42.69	45.93	34.06	42.61	76.01	26.83	55.64	49.25	46.47	STAT5B	signal transducer and activator of transcription 5B%2C transcript variant X1
5	508541	93.26	91.58	81.78	102.95	83.50	101.80	89.40	103.78	102.56	115.26	98.36	99.88	124.89	STAT3	signal transducer and activator of transcription 3
6	538961	10.22	8.55	9.42	9.27	9.08	8.02	14.13	7.54	9.37	9.48	6.36	13.40	9.46	ANKRD54	ankyrin repeat domain 54

$`GO:0007507 - heart development`

	GeneID	SLICK_2041	SLICK_2128	SLICK_2129	SLICK_2171	SLICK_2202	SLICK_2287	SLICK_466	WT_2058	WT_2068	WT_2089	WT_2116	WT_2226	WT_2273	gene	product
1	280939	16.94	14.04	20.79	11.52	16.07	14.08	15.52	13.27	21.43	11.18	29.90	18.78	16.79	TEK	TEK receptor tyrosine kinase
2	281134	1.86	3.77	6.49	2.30	5.05	6.06	4.63	1.66	4.81	2.00	2.56	3.78	3.95	ECE2	endothelin converting enzyme 2%2C transcript variant 2a
3	281193	7.28	10.21	6.89	10.15	5.37	7.58	6.03	10.50	7.01	12.63	8.67	5.40	6.48	GJA1	gap junction protein alpha 1
4	281270	2.70	1.46	1.59	4.59	1.34	2.90	3.14	3.79	1.18	4.91	3.98	5.02	3.18	HOPX	HOP homeobox%2C transcript variant X3
5	281733	5.07	3.13	3.26	8.76	4.79	7.57	6.95	9.17	4.72	5.98	5.80	4.56	9.08	CXADR	CXADR%2C Ig-like cell adhesion molecule%2C transcript variant X2
6	282321	12.91	13.52	6.38	13.61	9.82	12.60	13.43	16.22	11.38	13.84	15.60	8.61	11.73	PPP3R1	protein phosphatase 3 regulatory subunit B%2C alpha
7	282382	28.25	33.40	32.41	22.80	36.04	41.62	34.63	29.67	37.23	25.04	30.00	26.82	34.54	TGFBR1	transforming growth factor beta receptor 1%2C transcript variant X2
8	282682	96.48	118.55	104.65	58.00	113.50	115.97	96.90	82.98	118.92	63.92	101.75	145.76	110.85	GRK2	G protein-coupled receptor kinase 2%2C transcript variant X1
9	338068	13.62	12.41	7.58	17.67	10.74	12.82	10.20	19.40	10.65	21.30	14.80	9.09	13.11	ACVR1	activin A receptor type 1%2C transcript variant X1
10	506085	6.25	6.10	7.41	7.26	6.38	3.96	3.94	5.50	4.90	5.65	7.26	6.91	7.01	ZDHHC16	zinc finger DHHC-type containing 16%2C transcript variant X3
11	509656	13.93	20.62	19.93	13.90	17.87	16.26	16.60	15.71	19.94	14.95	19.50	22.53	19.89	PSKH1	protein serine kinase H1
12	512254	14.28	14.09	11.33	14.39	11.81	10.63	13.29	17.31	12.14	12.17	24.27	10.09	13.61	MEF2C	myocyte enhancer factor 2C%2C transcript variant X9
13	514997	8.78	9.24	15.24	7.71	8.89	10.06	7.33	5.61	11.24	7.86	5.13	8.99	8.30	ZNHIT1	zinc finger HIT-type containing 1
14	534041	11.77	12.38	18.45	10.11	11.45	14.24	16.55	15.36	19.48	10.40	12.63	19.27	19.33	MICAL2	microtubule associated monooxygenase%2C calponin and LIM domain containing 2%2C transcript variant X1
15	534069	34.58	37.30	23.76	38.46	29.75	29.27	27.73	48.82	27.82	39.82	32.13	33.18	33.27	TGFB2	transforming growth factor beta 2%2C transcript variant X3
16	538690	11.47	14.06	21.67	11.65	19.24	5.79	17.01	9.73	20.46	12.49	15.66	15.03	15.97	ID3	inhibitor of DNA binding 3%2C HLH protein
17	539430	17.64	22.10	19.65	15.66	27.24	15.05	19.45	19.32	18.72	20.67	10.97	14.14	8.76	OSR1	odd-skipped related transciption factor 1
18	539696	21.82	21.32	20.85	21.81	26.57	23.31	24.24	20.41	28.90	22.22	27.81	21.44	22.11	HEXIM1	hexamethylene bisacetamide inducible 1
19	615492	6.96	6.87	5.74	8.43	5.27	5.14	8.09	4.82	9.53	8.39	4.00	7.62	5.34	SMG9	SMG9%2C nonsense mediated mRNA decay factor
20	786335	18.30	16.31	14.01	20.27	16.04	21.45	12.42	22.04	17.71	21.72	15.49	15.65	22.29	GYS1	glycogen synthase 1

$`GO:0009409 - response to cold`

	GeneID	SLICK_2041	SLICK_2128	SLICK_2129	SLICK_2171	SLICK_2202	SLICK_2287	SLICK_466	WT_2058	WT_2068	WT_2089	WT_2116	WT_2226	WT_2273	gene	product
1	281605	2.75	2.14	1.44	2.01	3.18	1.46	4.15	2.56	2.82	2.21	3.29	3.67	5.86	ADRB2	adrenoceptor beta 2
2	281827	0.89	1.44	1.21	0.77	0.41	0.57	0.82	1.34	1.12	1.91	0.78	0.47	1.00	HSPA2	heat shock protein family A (Hsp70) member 2
3	281832	238.29	220.47	234.09	303.29	249.61	302.17	218.58	360.65	258.21	380.33	307.61	202.29	234.78	HSP90AA1	heat shock protein 90 alpha family class A member 1
4	282130	45.14	47.64	56.03	49.58	55.77	40.32	46.71	36.15	64.45	47.91	49.00	69.18	53.79	ACADVL	acyl-CoA dehydrogenase very long chain%2C transcript variant X1
5	505968	32.49	42.72	26.03	18.63	25.97	25.45	21.49	31.01	25.94	23.17	46.37	34.32	25.33	ACADM	acyl-CoA dehydrogenase medium chain%2C transcript variant X1
6	523885	8.66	11.95	12.11	7.24	9.13	13.89	12.79	10.67	12.15	8.50	9.96	11.69	17.61	TMEM135	transmembrane protein 135%2C transcript variant X2

$`GO:0001942 - hair follicle development`

	GeneID	SLICK_2041	SLICK_2128	SLICK_2129	SLICK_2171	SLICK_2202	SLICK_2287	SLICK_466	WT_2058	WT_2068	WT_2089	WT_2116	WT_2226	WT_2273	gene	product
1	281867	0.19	0.21	0.24	0.40	0.32	0.12	0.08	0.48	1.34	1.31	0.42	0.76	0.44	INHBA	inhibin subunit beta A%2C transcript variant X2
2	505423	50.79	35.89	24.35	44.42	29.23	38.92	36.78	61.48	35.73	45.27	61.38	31.31	37.88	LGR4	leucine rich repeat containing G protein-coupled receptor 4
3	534069	34.58	37.30	23.76	38.46	29.75	29.27	27.73	48.82	27.82	39.82	32.13	33.18	33.27	TGFB2	transforming growth factor beta 2%2C transcript variant X3
4	535814	8.58	8.86	7.08	8.28	8.08	9.00	8.56	11.05	6.74	6.76	7.93	5.55	7.27	ZDHHC21	zinc finger DHHC-type containing 21%2C transcript variant X1
5	616694	18.67	35.43	59.20	13.69	35.40	28.02	28.15	24.75	44.91	15.86	25.00	48.52	37.63	NSDHL	NAD(P) dependent steroid dehydrogenase-like%2C transcript variant X1

$`GO:0097009 - energy homeostasis`

	GeneID	SLICK_2041	SLICK_2128	SLICK_2129	SLICK_2171	SLICK_2202	SLICK_2287	SLICK_466	WT_2058	WT_2068	WT_2089	WT_2116	WT_2226	WT_2273	gene	product
1	281052	1147.38	1848.11	2525.53	393.44	2177.52	2320.86	1873.41	1606.57	2062.02	603.62	3116.64	2527.42	2110.53	CD36	CD36 molecule%2C transcript variant 1
2	281370	312.28	445.53	347.16	276.91	365.57	392.35	359.19	309.06	341.74	260.43	434.19	471.35	395.46	UBB	ubiquitin B
3	282306	18.98	19.75	14.58	16.91	18.76	16.11	19.08	36.41	20.04	21.77	26.90	9.71	18.77	PIK3CA	phosphatidylinositol-4%2C5-bisphosphate 3-kinase catalytic subunit alpha%2C transcript variant X1
4	338446	0.47	12.26	6.83	0.77	16.38	3.67	3.49	8.19	12.32	0.35	1.55	13.73	1.79	PPARGC1A	PPARG coactivator 1 alpha
5	506223	2.94	2.62	2.20	2.42	1.52	3.00	1.90	2.37	2.20	2.73	3.01	1.40	4.17	OMA1	OMA1 zinc metallopeptidase%2C transcript variant X5
6	508541	93.26	91.58	81.78	102.95	83.50	101.80	89.40	103.78	102.56	115.26	98.36	99.88	124.89	STAT3	signal transducer and activator of transcription 3
7	513852	0.73	0.00	0.69	0.54	0.45	0.65	1.07	0.47	0.89	1.08	1.05	1.37	1.17	PM20D1	peptidase M20 domain containing 1
8	526652	15.14	12.79	10.45	17.93	13.73	14.27	14.97	18.36	11.11	15.96	16.19	12.94	13.07	FLCN	folliculin%2C transcript variant X1

$`GO:0051260 - protein homooligomerization`

	GeneID	SLICK_2041	SLICK_2128	SLICK_2129	SLICK_2171	SLICK_2202	SLICK_2287	SLICK_466	WT_2058	WT_2068	WT_2089	WT_2116	WT_2226	WT_2273	gene	product
1	100124502	4.22	3.19	3.50	4.47	5.48	3.89	4.59	5.38	3.73	3.62	3.09	2.93	2.70	KCTD7	potassium channel tetramerization domain containing 7
2	100125308	15.63	15.99	20.72	9.89	17.84	16.69	15.11	11.64	13.71	9.22	10.18	20.56	20.99	KCTD5	potassium channel tetramerization domain containing 5%2C transcript variant X1
3	281088	0.28	0.00	0.34	1.05	0.11	0.69	0.25	0.30	0.35	2.61	0.00	0.47	0.39	COMP	cartilage oligomeric matrix protein
4	281426	0.13	0.00	0.11	2.76	0.04	0.49	0.04	0.80	0.42	0.35	0.00	0.00	0.44	PRND	prion like protein doppel
5	281427	55.69	72.44	50.61	38.97	68.48	68.99	61.35	97.43	52.71	69.51	64.00	41.92	64.76	PRNP	prion protein%2C transcript variant 1
6	281465	4.51	3.14	3.61	5.10	4.54	4.10	5.33	1.94	3.09	3.95	3.48	3.60	3.13	ROM1	retinal outer segment membrane protein 1%2C transcript variant X1
7	281528	16.30	20.89	34.81	24.19	24.65	15.61	15.33	22.96	32.06	22.24	22.46	35.60	19.73	TGM2	transglutaminase 2
8	281880	1.01	0.00	0.11	1.37	0.30	0.61	0.33	0.72	0.54	1.34	0.83	0.00	0.17	KCNA4	potassium voltage-gated channel subfamily A member 4
9	282573	14.11	14.90	28.08	11.81	18.83	20.86	15.45	14.48	14.49	12.91	15.57	42.96	14.94	KCNMA1	potassium calcium-activated channel subfamily M alpha 1%2C transcript variant X8
10	282846	1.04	2.31	1.26	3.10	2.07	1.24	1.13	1.73	0.66	2.51	0.58	3.02	1.62	PYCARD	PYD and CARD domain containing
11	317724	63.33	63.55	70.56	37.98	71.49	92.57	71.51	74.74	73.12	31.26	95.55	75.67	66.59	USO1	USO1 vesicle transport factor%2C transcript variant X1
12	507899	26.41	24.50	29.39	23.86	26.68	25.08	19.47	28.84	23.62	26.38	28.02	22.87	26.28	MICU1	mitochondrial calcium uptake 1
13	507911	1.68	2.35	2.36	2.50	1.23	1.17	1.03	1.04	1.54	2.04	1.01	1.37	1.34	KCTD13	potassium channel tetramerization domain containing 13
14	508343	4.62	6.25	4.73	3.32	7.64	5.42	4.70	3.41	6.28	4.76	5.51	4.61	3.16	JMJD6	arginine demethylase and lysine hydroxylase
15	510679	17.90	22.06	19.37	23.63	21.54	16.50	22.30	11.86	17.23	17.69	17.39	22.01	17.24	ALAD	aminolevulinate dehydratase
16	511908	15.93	25.69	19.13	24.59	19.81	18.61	14.80	18.54	19.92	24.43	13.35	17.14	21.06	EHD1	EH domain containing 1
17	512402	18.13	20.70	23.04	24.02	19.92	16.57	17.32	15.26	17.36	20.75	13.41	17.42	18.00	SHKBP1	SH3KBP1 binding protein 1%2C transcript variant X2
18	512480	0.65	0.80	0.56	0.75	0.69	0.74	0.72	0.51	0.21	0.73	0.20	0.17	0.86	NLRC4	NLR family CARD domain containing 4%2C transcript variant X2
19	512578	5.75	2.71	3.42	4.34	4.83	2.87	6.62	3.37	6.29	4.69	2.73	2.21	1.86	KCTD15	potassium channel tetramerization domain containing 15%2C transcript variant X1
20	513324	40.08	62.41	84.54	30.66	62.29	50.53	54.70	33.16	67.09	31.07	43.23	82.02	52.01	LETM1	leucine zipper and EF-hand containing transmembrane protein 1%2C transcript variant X1
21	518384	3.40	4.19	1.60	3.28	2.12	2.74	2.35	3.37	1.60	2.66	2.12	2.59	1.81	KCND1	potassium voltage-gated channel subfamily D member 1%2C transcript variant X2
22	520173	17.05	19.64	19.47	14.58	19.26	17.41	15.11	13.18	18.26	10.22	14.22	22.61	19.67	TMEM120A	transmembrane protein 120A
23	521442	10.72	12.68	19.76	8.79	15.49	16.42	14.38	11.79	13.19	9.98	11.52	16.75	12.43	SPAST	spastin%2C transcript variant X2
24	521998	10.39	12.90	13.57	7.72	14.67	12.71	12.72	9.29	12.36	9.28	10.72	18.29	11.38	ZNF746	zinc finger protein 746%2C transcript variant X3
25	535424	0.26	0.14	0.06	0.39	0.56	0.73	0.12	0.83	0.04	0.26	0.49	0.31	0.00	ATL1	atlastin GTPase 1%2C transcript variant X2
26	539167	2.85	3.08	2.99	5.19	2.20	3.16	2.83	3.56	3.20	4.62	3.26	2.15	3.04	KCTD11	potassium channel tetramerization domain containing 11
27	613694	4.89	4.27	1.81	5.53	1.94	3.31	5.05	5.82	3.13	3.47	5.60	3.00	4.36	KCTD18	potassium channel tetramerization domain containing 18%2C transcript variant X1
28	783855	9.24	8.29	10.66	5.68	8.04	8.83	10.17	9.28	7.10	5.45	10.23	8.37	8.83	TIFA	TRAF interacting protein with forkhead associated domain%2C transcript variant X1
29	784587	11.86	13.63	10.47	22.04	11.13	12.58	13.57	22.65	8.88	23.00	16.32	7.72	12.33	KCTD1	potassium channel tetramerization domain containing 1%2C transcript variant X3

$`GO:0140668 - positive regulation of oxytocin production`

	GeneID	SLICK_2041	SLICK_2128	SLICK_2129	SLICK_2171	SLICK_2202	SLICK_2287	SLICK_466	WT_2058	WT_2068	WT_2089	WT_2116	WT_2226	WT_2273	gene	product
1	281416	4.58	5.26	3.22	2.43	4.83	3.59	3.51	1.86	2.83	0.85	1.18	4.44	2.46	POMC	proopiomelanocortin%2C transcript variant X1

$`GO:0060397 - growth hormone receptor signaling pathway via JAK-STAT`

	GeneID	SLICK_2041	SLICK_2128	SLICK_2129	SLICK_2171	SLICK_2202	SLICK_2287	SLICK_466	WT_2058	WT_2068	WT_2089	WT_2116	WT_2226	WT_2273	gene	product
1	282376	48.07	59.84	50.23	25.17	42.69	45.93	34.06	42.61	76.01	26.83	55.64	49.25	46.47	STAT5B	signal transducer and activator of transcription 5B%2C transcript variant X1
2	508541	93.26	91.58	81.78	102.95	83.50	101.80	89.40	103.78	102.56	115.26	98.36	99.88	124.89	STAT3	signal transducer and activator of transcription 3

$`GO:0006816 - calcium ion transport`

	GeneID	SLICK_2041	SLICK_2128	SLICK_2129	SLICK_2171	SLICK_2202	SLICK_2287	SLICK_466	WT_2058	WT_2068	WT_2089	WT_2116	WT_2226	WT_2273	gene	product
1	281040	134.56	171.76	163.81	94.97	166.61	100.67	109.04	165.08	171.85	135.34	243.58	165.44	145.98	CAV1	caveolin 1%2C transcript variant X1
2	281879	13.99	7.12	7.20	13.15	5.92	4.15	10.65	8.38	10.59	15.81	8.20	6.84	8.22	ITPR3	inositol 1%2C4%2C5-trisphosphate receptor type 3%2C transcript variant X1
3	282178	0.63	0.45	0.40	0.21	0.22	0.00	1.03	0.52	0.27	0.39	0.26	1.17	1.25	CHRNA7	cholinergic receptor%2C nicotinic%2C alpha 7
4	282196	22.30	12.54	8.41	35.64	12.81	10.16	14.50	16.71	10.70	21.22	11.06	12.95	14.31	CORO1A	coronin 1A%2C transcript variant X2
5	282325	7.09	5.41	3.39	12.47	3.55	9.06	7.05	10.07	2.26	6.57	6.20	3.93	6.04	PRKCB	protein kinase C beta%2C transcript variant X1
6	327663	59.22	48.03	55.32	41.95	55.45	50.47	51.17	57.10	56.41	44.78	57.22	50.56	63.71	ATP2C1	ATPase secretory pathway Ca2+ transporting 1%2C transcript variant X3
7	327685	22.98	28.36	28.64	46.71	25.36	21.45	25.80	31.97	32.45	43.43	21.94	22.98	23.74	ANXA6	annexin A6%2C transcript variant X1
8	337925	13.37	15.87	24.54	8.52	30.12	28.53	23.02	18.03	27.59	12.58	27.47	27.74	26.95	SLC8A1	solute carrier family 8 member A1%2C transcript variant X3
9	515461	121.43	129.84	163.84	85.40	145.33	177.80	113.42	100.39	131.44	89.99	74.01	110.99	113.95	SARAF	store-operated calcium entry associated regulatory factor
10	540892	16.55	17.64	18.58	15.98	18.78	26.54	18.44	20.31	20.44	15.35	22.42	24.04	26.10	DNM1L	dynamin 1 like%2C transcript variant X7

$`GO:0050830 - defense response to Gram-positive bacterium`

	GeneID	SLICK_2041	SLICK_2128	SLICK_2129	SLICK_2171	SLICK_2202	SLICK_2287	SLICK_466	WT_2058	WT_2068	WT_2089	WT_2116	WT_2226	WT_2273	gene	product
1	281249	1.21	1.67	0.24	4.56	0.08	1.88	1.14	1.99	0.23	1.49	0.79	0.70	0.22	IL18	interleukin 18%2C transcript variant X1
2	338039	17.13	16.35	13.41	21.30	12.37	21.55	13.42	15.23	11.30	13.38	14.60	10.37	20.35	CASP4	caspase 4%2C apoptosis-related cysteine peptidase%2C transcript variant X1
3	444881	33.75	27.58	56.08	21.93	49.11	46.53	27.28	21.91	35.02	21.10	32.08	37.67	48.47	MYD88	myeloid differentiation primary response 88%2C transcript variant X2
4	506206	1.06	0.28	0.26	1.65	0.69	0.88	0.83	1.33	0.32	0.88	0.71	0.00	0.00	MR1	major histocompatibility complex%2C class I-related%2C transcript variant X3
5	506509	22.63	22.54	27.20	19.75	18.56	32.85	26.56	23.95	30.55	18.16	28.05	22.37	23.35	SEH1L	SEH1 like nucleoporin%2C transcript variant X1
6	508990	12.30	17.91	22.01	26.94	16.62	10.93	19.07	9.63	22.20	19.54	9.49	14.59	14.10	RARRES2	retinoic acid receptor responder 2%2C transcript variant X2
7	534407	14.05	12.36	13.89	12.36	14.44	16.99	15.96	16.17	14.67	10.10	17.19	13.67	14.85	RIPK2	receptor interacting serine/threonine kinase 2%2C transcript variant X2
8	539997	14.15	14.82	9.13	38.25	17.06	14.68	8.85	14.57	11.54	24.29	10.50	14.56	18.21	MPEG1	macrophage expressed 1
9	540444	16.31	9.95	9.77	21.00	13.12	10.93	13.01	13.09	11.32	19.31	17.39	12.11	14.85	HMGB2	high mobility group box 2
10	618530	10.73	9.80	15.03	9.68	15.02	9.21	9.60	5.92	19.50	10.74	7.70	14.11	10.02	ROMO1	reactive oxygen species modulator 1%2C transcript variant X1
11	767908	155.49	150.63	143.32	142.40	168.43	135.96	175.34	112.17	164.50	124.46	101.98	144.32	135.06	RPL39	ribosomal protein L39

$`GO:0003085 - negative regulation of systemic arterial blood pressure`

	GeneID	SLICK_2041	SLICK_2128	SLICK_2129	SLICK_2171	SLICK_2202	SLICK_2287	SLICK_466	WT_2058	WT_2068	WT_2089	WT_2116	WT_2226	WT_2273	gene	product
1	527572	3.64	1.48	4.37	1.28	4.45	1.66	0.53	0.15	0.69	1.68	0.34	1.02	1.86	GPR37L1	G protein-coupled receptor 37 like 1
2	532120	23.25	29.76	17.70	16.74	21.53	25.34	24.55	23.78	23.18	14.61	19.48	23.65	26.04	BBS4	Bardet-Biedl syndrome 4%2C transcript variant X4

$`GO:0042310 - vasoconstriction`

	GeneID	SLICK_2041	SLICK_2128	SLICK_2129	SLICK_2171	SLICK_2202	SLICK_2287	SLICK_466	WT_2058	WT_2068	WT_2089	WT_2116	WT_2226	WT_2273	gene	product
1	281040	134.56	171.76	163.81	94.97	166.61	100.67	109.04	165.08	171.85	135.34	243.58	165.44	145.98	CAV1	caveolin 1%2C transcript variant X1
2	527114	0.54	0.70	1.78	0.45	0.75	3.19	3.13	0.92	5.91	1.15	0.89	2.93	1.69	AGT	angiotensinogen

$`GO:0051302 - regulation of cell division`

	GeneID	SLICK_2041	SLICK_2128	SLICK_2129	SLICK_2171	SLICK_2202	SLICK_2287	SLICK_466	WT_2058	WT_2068	WT_2089	WT_2116	WT_2226	WT_2273	gene	product
1	327682	579.77	592.48	493.33	547.93	527.13	665.62	540.97	639.48	438.55	618.07	652.12	389.17	689.11	RACK1	receptor for activated C kinase 1
2	505401	4.89	4.92	6.94	8.65	5.41	6.07	7.13	7.11	6.01	7.92	6.32	7.95	6.10	TADA2A	transcriptional adaptor 2A%2C transcript variant X3
3	507401	1.76	0.85	1.55	1.27	0.79	1.72	0.38	0.57	0.12	0.56	0.29	1.12	0.85	OOEP	oocyte expressed protein
4	510124	6.59	11.27	8.64	4.12	10.79	6.44	7.65	5.62	8.57	4.87	10.32	7.86	5.23	EFHC1	EF-hand domain containing 1
5	510141	17.39	16.40	23.40	17.63	17.32	14.15	14.69	11.56	18.94	15.08	12.81	18.76	20.63	CIB1	calcium and integrin binding 1
6	510184	21.50	18.12	18.12	29.04	22.47	21.98	23.95	20.03	18.79	26.09	17.26	16.10	16.25	TADA3	transcriptional adaptor 3%2C transcript variant X1

Figure

counts.per.million <- cpm(y) %>%
  as.data.frame %>%
  rownames_to_column("GeneID")

final.go.terms <- data.frame(GOID=sprintf("GO:%07d", c(60397, 9409, 1942, 42310,
                                       140668, 51260, 50830)))
plot.names <- go.fry %>%
  rownames_to_column('GOID') %>%
  left_join(go.names, by='GOID') %>%
  right_join(final.go.terms) %>%
  split(.$GOID)

go.data <- lapply(plot.names,
                  function(x) filter(counts.per.million,
                                     GeneID %in% go.sets[[x$GOID]]))

go.plots <- mapply(function(data, name) {
  gather(data, -GeneID, key='sample', value='cpm') %>%
    left_join(desc, by="GeneID")  %>%
    mutate(group=sub("_.*", "", sample),
           product = sub(', transcript variant X[1-9]', '', URLdecode(product))) %>%
    ggplot(aes(sample, product, alpha=cpm, fill=group)) +
    geom_tile(color='grey') +
    geom_text(aes(label=sprintf("%0.1f", cpm)), alpha=1)+
    scale_fill_manual(values=c(WT="#cc546e", SLICK="#4db598"), labels=c(WT="Wild-Type", SLICK="SLICK")) +
    scale_alpha_continuous() +
    ggtitle(sprintf("%s: %s (%s)", name$GOID, name$TERM, name$Direction)) +
    theme_minimal() +
    theme(panel.grid = element_blank(),
          legend.position = 'none',
          axis.title =element_blank(),
          axis.text.x = element_blank())
}, go.data, plot.names, SIMPLIFY = F)

library(cowplot)

short.plots <- go.plots[names(go.plots) != "GO:0051260"]
rel.heights <- sapply(go.data[names(go.data) != "GO:0051260"], nrow) + 1.5

plot_grid(go.plots[["GO:0051260"]] + theme(axis.text.x=element_text(angle=45, hjust=1, vjust=1)),
          plot_grid(plotlist=short.plots, ncol=1, rel_heights = rel.heights , align="v" ),
          ncol=2)

Session Info

sessionInfo()

Meta - Reduced

library(Hmisc)
library(tidyverse)
library(cowplot)

metabolites <-
  read.delim("metabolomics.reduced.tbl", nrows = 1) %>% select(-c(1:8)) %>% c %>% unlist

data <- read.delim("metabolomics.reduced.tbl", skip=1) %>%
  rename(setNames(make.names(metabolites), names(metabolites)),
         Sample="Sample..")

meta.data <- data[,c(1:3)]
phys <- data %>%
  select(1, c(4:8)) %>%
  column_to_rownames("Sample") %>%
  mutate_all(as.numeric)
metab <- data %>%
  select("Sample", names(metabolites)) %>%
  column_to_rownames("Sample") %>%
  mutate_all(as.numeric)

cor.data <- rcorr(as.matrix(metab), as.matrix(phys), type="spearman") %>%
  lapply('[', names(metabolites), names(phys)) %>%
  lapply(as.data.frame) %>%
  lapply(rownames_to_column, "ID") %>%
  bind_rows(.id='Type') %>%
  gather( "Variable", "val", -ID, -Type ) %>%
  spread(Type, val) %>%
  mutate(color = fct_recode(cut(P, breaks=c(-Inf, 0.01, 0.05 , Inf))),
         Variable = factor(Variable, levels=c("MY", "VT", "RI", "Diameter", "Total_MBF")),
         ID = factor(ID, levels=names(sort(metabolites))),
         set = (as.numeric(ID)+1) %% 2)

levels(cor.data$color) <- c("<=0.01", "0.01<x<=0.05", ">0.05")

metabolite.labels <- setNames(sprintf("%s\n%s", names(sort(metabolites)), sort(metabolites)), names(sort(metabolites)))
plot <- ggplot(cor.data, aes(Variable, ID,
                     size=abs(r), alpha=abs(r), fill=r,
                     color=color)) +
  geom_point(shape=21, stroke=2) +
  scale_x_discrete(position = 'top') +
  scale_y_discrete(limits=rev, labels=metabolite.labels) +
  scale_alpha_continuous(guide='none', range=c(1,1), limits=c(0,1)) +
  scale_size_continuous(guide='none', range=c(1,15), limits=c(0,1)) +
  scale_fill_distiller( palette = 'RdBu', direction = 1,
                        guide = guide_colorbar(barwidth=25), name='Correlation', limits=c(-1,1)) +
  scale_color_manual(values=c("#D4A10D", "#848284", 'white'), name="Significance") +
  facet_wrap(set~., scales="free", ncol=2) +
  theme_minimal() +
  theme(axis.title = element_blank(),
        panel.grid = element_blank(),
        strip.background = element_blank(),
        strip.text.x = element_blank(),
        legend.position = 'none',
        legend.box="vertical")

plot_grid(plot, get_legend(plot + theme(legend.position = 'bottom')), nrow = 2, rel_heights = c(1, 0.1))

options(asciiType = "org")
mutate(cor.data, value = sprintf("%.03f (%.03f)", r, P)) %>%
  select(ID, Variable, value) %>%
  spread("Variable", "value") %>%
  mutate(ID = factor(ID, levels=names(sort(metabolites)))) %>%
  ascii::ascii()

	ID	MY	VT	RI	Diameter	Total_MBF
1	HMDB0244503	-0.617 (0.019)	0.538 (0.047)	0.503 (0.067)	-0.253 (0.383)	-0.433 (0.122)
2	HMDB0003175	0.582 (0.029)	-0.358 (0.208)	-0.090 (0.759)	0.451 (0.106)	0.262 (0.366)
3	HMDB0006245	-0.337 (0.238)	0.174 (0.553)	0.604 (0.022)	-0.613 (0.020)	-0.574 (0.032)
4	HMDB0244869	0.710 (0.004)	0.020 (0.946)	-0.477 (0.085)	0.380 (0.180)	0.530 (0.051)
5	HMDB0245244	-0.412 (0.143)	0.266 (0.358)	0.705 (0.005)	-0.495 (0.072)	-0.569 (0.034)
6	HMDB0037544	-0.710 (0.004)	0.332 (0.246)	0.319 (0.267)	-0.622 (0.018)	-0.547 (0.043)
7	HMDB0000528	0.822 (0.000)	-0.196 (0.503)	-0.204 (0.483)	0.125 (0.670)	0.332 (0.246)
8	HMDB0062738	0.406 (0.150)	-0.138 (0.637)	-0.565 (0.035)	0.516 (0.059)	0.785 (0.001)
9	HMDB0040551	0.157 (0.593)	-0.235 (0.418)	-0.644 (0.013)	0.569 (0.034)	0.543 (0.045)
10	HMDB0304536	-0.348 (0.222)	0.121 (0.681)	0.407 (0.149)	-0.516 (0.059)	-0.640 (0.014)
11	HMDB0002212	-0.774 (0.001)	0.059 (0.840)	0.415 (0.140)	-0.420 (0.135)	-0.451 (0.106)
12	HMDB0000063	-0.042 (0.887)	0.275 (0.342)	0.015 (0.958)	-0.398 (0.159)	-0.538 (0.047)
13	HMDB0242441	0.077 (0.793)	-0.732 (0.003)	-0.380 (0.180)	0.385 (0.175)	0.182 (0.533)
14	HMDB0001051	-0.329 (0.251)	0.284 (0.326)	0.697 (0.006)	-0.411 (0.144)	-0.569 (0.034)
15	HMDB0000131	-0.013 (0.964)	0.666 (0.009)	0.323 (0.260)	-0.503 (0.067)	-0.266 (0.358)
16	HMDB0002345	-0.664 (0.010)	0.459 (0.098)	0.495 (0.072)	-0.354 (0.215)	-0.525 (0.054)
17	HMDB0010381	0.545 (0.044)	-0.160 (0.584)	-0.099 (0.737)	0.301 (0.296)	0.442 (0.114)
18	HMDB0010390	0.404 (0.152)	-0.688 (0.007)	-0.530 (0.051)	0.376 (0.185)	0.349 (0.221)
19	HMDB0000156	-0.236 (0.417)	0.367 (0.197)	0.675 (0.008)	-0.565 (0.035)	-0.569 (0.034)
20	HMDB0255110	-0.659 (0.010)	0.407 (0.149)	0.600 (0.023)	-0.468 (0.091)	-0.508 (0.064)
21	HMDB0062342	0.181 (0.536)	-0.710 (0.004)	-0.385 (0.175)	0.178 (0.543)	0.178 (0.543)
22	HMDB0002088	-0.670 (0.009)	0.420 (0.135)	0.516 (0.059)	-0.495 (0.072)	-0.530 (0.051)
23	HMDB0013648	-0.461 (0.097)	0.341 (0.233)	0.459 (0.098)	-0.675 (0.008)	-0.670 (0.009)
24	HMDB0002100	-0.723 (0.003)	0.376 (0.185)	0.376 (0.185)	-0.477 (0.085)	-0.451 (0.106)
25	HMDB0029006	0.549 (0.042)	-0.209 (0.474)	-0.468 (0.091)	0.578 (0.030)	0.530 (0.051)
26	HMDB0002752	-0.364 (0.201)	0.037 (0.899)	0.200 (0.493)	-0.626 (0.017)	-0.596 (0.025)
27	HMDB0001403	0.598 (0.024)	-0.398 (0.159)	-0.182 (0.533)	0.301 (0.296)	0.262 (0.366)
28	HMDB0001220	0.426 (0.129)	-0.327 (0.253)	0.169 (0.563)	0.600 (0.023)	0.222 (0.446)
29	HMDB0002664	0.659 (0.010)	-0.108 (0.714)	-0.073 (0.805)	0.481 (0.081)	0.415 (0.140)
30	HMDB0030569	0.289 (0.317)	-0.495 (0.072)	-0.618 (0.019)	0.112 (0.703)	0.244 (0.401)
31	HMDB0003252	0.602 (0.023)	-0.310 (0.281)	-0.169 (0.563)	0.297 (0.303)	0.385 (0.175)
32	HMDB0014447	-0.234 (0.421)	0.402 (0.154)	0.679 (0.008)	-0.596 (0.025)	-0.648 (0.012)
33	HMDB0000929	-0.157 (0.593)	0.473 (0.088)	0.578 (0.030)	-0.407 (0.149)	-0.275 (0.342)
34	HMDB0029090	0.452 (0.105)	-0.182 (0.533)	-0.301 (0.296)	0.640 (0.014)	0.398 (0.159)

Name		Name	Last commit message	Last commit date
Latest commit History 8 Commits
0-ref		0-ref
1-raw		1-raw
2-salmon		2-salmon
apps		apps
.gitignore		.gitignore
BCV.png		BCV.png
BCV.rm.png		BCV.rm.png
MDS.png		MDS.png
MDS.rm.png		MDS.rm.png
bta00340.pathview.png		bta00340.pathview.png
bta00515.pathview.png		bta00515.pathview.png
bta00562.pathview.png		bta00562.pathview.png
bta00590.pathview.png		bta00590.pathview.png
bta00591.pathview.png		bta00591.pathview.png
bta04012.pathview.png		bta04012.pathview.png
bta04070.pathview.png		bta04070.pathview.png
bta04530.pathview.png		bta04530.pathview.png
go-genes.xlsx		go-genes.xlsx
heatmap.go.png		heatmap.go.png
mean-diff.lrt.png		mean-diff.lrt.png
mean-diff.lrt.rm.png		mean-diff.lrt.rm.png
mean-diff.png		mean-diff.png
mean-diff.rm.png		mean-diff.rm.png
meta.corr.png		meta.corr.png
meta.data.tbl		meta.data.tbl
meta.reduced.corr.png		meta.reduced.corr.png
metabolomics.reduced.tbl		metabolomics.reduced.tbl
readme.org		readme.org
report.org		report.org
venn-MDS.png		venn-MDS.png
venn-MDS.rm.png		venn-MDS.rm.png
venn.png		venn.png
venn.rm.png		venn.rm.png
volcano.lrt.png		volcano.lrt.png
volcano.lrt.rm.png		volcano.lrt.rm.png

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Repository files navigation

SLICK Holstein mammary gland RNA-Seq methods

Bos taurus reference

Raw data

Salmon

EdgeR

Load data

Normalize counts, Filter Genes, and Estimate Dispersions

DEG LRT test

KEGG

GO

Session Info

EdgeR w/o Outlier

Load data

Normalize counts, Filter Genes, and Estimate Dispersions

DEG LRT test

KEGG

GO

GO terms of interest gene expression

Figure

Session Info

Meta - Reduced

About

IGBB/slick.rna-seq

Folders and files

Latest commit

History

Repository files navigation

SLICK Holstein mammary gland RNA-Seq methods

Bos taurus reference

Raw data

Salmon

EdgeR

Load data

Normalize counts, Filter Genes, and Estimate Dispersions

DEG LRT test

KEGG

GO

Session Info

EdgeR w/o Outlier

Load data

Normalize counts, Filter Genes, and Estimate Dispersions

DEG LRT test

KEGG

GO

GO terms of interest gene expression

Figure

Session Info

Meta - Reduced

About

Resources

Stars

Watchers

Forks