Raphaël C1q data analysis

This page described the code and the result of C1q data shared by Inge.

Identification of differential acceessibilities

First, the region x sample count matrix was loaded in R and removed non-acceessible regions in most of the samples

inge = read.table("ChIP-seq raw technical reps n4 RPMI C1q.txt")
library(edgeR)
cpms = cpm(inge)
keep = rowSums(cpms > 1) >= 2
inge = inge[keep,]

Because this is a multifactor design and donor variability should be removed by GLM. Therefore, we make a design matrix and donor-variability is removed.

coldata =data.frame(row.names = colnames(inge))
coldata$donor = c(rep("donor1",2),rep("donor2",2),rep("donor3",2),rep("donor4",2))
coldata$condition = rep(c("RPMI","C1q"),4)
design = model.matrix( ~ donor + condition, coldata)
design

       (Intercept) donordonor2 donordonor3 donordonor4 conditionRPMI
RPMIn1           1           0           0           0             1
C1qn1            1           0           0           0             0
RPMIn2           1           1           0           0             1
C1qn2            1           1           0           0             0
RPMIn3           1           0           1           0             1
C1qn3            1           0           1           0             0
RPMIn4           1           0           0           1             1
C1qn4            1           0           0           1             0
attr(,"assign")
[1] 0 1 1 1 2
attr(,"contrasts")
attr(,"contrasts")$donor
[1] "contr.treatment"

attr(,"contrasts")$condition
[1] "contr.treatment"

d = DGEList(counts = inge, group = rep(c("RPMI","C1q"),4),sample = coldata)
d = calcNormFactors(d)
plotMDS(d, labels = rownames(coldata),col = c("darkgreen","blue")[factor(coldata$condition)])

d2 = estimateGLMTrendedDisp(d, design)
d2 = estimateGLMTagwiseDisp(d2, design)
f = glmFit(d2, design)
de = glmLRT(f, coef = 5)
tt = topTags(de, n = nrow(d))
head(tt$table)

Keep in mind that C1q is 0 and RPMI is 1 in this design matrix and this means that logFC > 0 is downregulated (i.e. this program calculats RMPI/C1q, not C1q/RPMI)

A good thing is that C1q and RPMI samples are likely clustered seperately in the MDS plot. But, the result suggests that the most of regions are not FDR<0.05. Now we checked the number of differential accessibilities in various threshold.

library(ggplot2)
das = tt$table
rownames(das) = paste0("chr",rownames(das))

#FDR < 0.05
tmp1 = nrow(subset(das,FDR < 0.05 & logFC < 0))
tmp2 = nrow(subset(das,FDR < 0.05 & logFC > 0))

#Pvalue < 0.01
tmp3 = nrow(subset(das,PValue < 0.01 & logFC < 0))
tmp4 = nrow(subset(das,PValue < 0.01 & logFC > 0))

#Pvalue < 0.05
tmp5 = nrow(subset(das,PValue < 0.05 & logFC < 0))
tmp6 = nrow(subset(das,PValue < 0.05 & logFC > 0))

df = data.frame(threshold=rep(c("FDR0.05","Pval0.01","Pval0.05"),each=2),Up_Down=rep(c("Up","Down"),3),Number=c(tmp1,-tmp2,tmp3,-tmp4,tmp5,-tmp6)) 
df$Up_Down = factor(df$Up_Down,levels = c("Up","Down"))

p = ggplot(df,aes(x= threshold,y=Number,fill=Up_Down)) + geom_bar(stat = "identity")
p + theme_minimal() + coord_flip()

As FDR< 0.05 threshold gives us almost no regions, we set P val <0.01 as a threshold. We may not be able to say “significantly” differential accessibilities, but we will see a trend of differential accessibilities induced by C1q.

Peak annotation

Annotation of all peaks

Now we identified differential accessibility between RPMI and C1q samples. Next step is to annotate peaks. To annotate, I use ChIPseeker. The coordinate notion in the original data is without “chr” and some software does not accept this format. So, I add “chr” in the coordinate before starting annotation. First, we look at all differential accessibilities (Upregulated and downregulated) and then check up- and down-regulated regions indivisually.

Here I show the peaks’ genomic context and some enriched ontologies in the text format and dotplot.

library(GenomicRanges)
library(ChIPseeker)
library(TxDb.Hsapiens.UCSC.hg38.knownGene)
library(org.Hs.eg.db)
library(clusterProfiler)
das = subset(tt$table,  PValue  < 0.01)
rownames(das) =paste0("chr",rownames(das))
gr = GRanges(rownames(das))
gr_df = as.data.frame(gr)
peakAnno <- annotatePeak(gr, tssRegion=c(-3000, 3000),TxDb=TxDb.Hsapiens.UCSC.hg38.knownGene, annoDb="org.Hs.eg.db")

>> preparing features information...         2024-01-29 14:11:46 
>> identifying nearest features...       2024-01-29 14:11:46 
>> calculating distance from peak to TSS...  2024-01-29 14:11:47 
>> assigning genomic annotation...       2024-01-29 14:11:47 
>> adding gene annotation...             2024-01-29 14:11:51 

'select()' returned 1:many mapping between keys and columns

>> assigning chromosome lengths          2024-01-29 14:11:51 
>> done...                   2024-01-29 14:11:52 

plotAnnoPie(peakAnno)

gene <- seq2gene(gr, tssRegion = c(-1000, 1000), flankDistance = 3000, TxDb=TxDb.Hsapiens.UCSC.hg38.knownGene)
pathway2 = enrichGO(gene,OrgDb = org.Hs.eg.db,ont="BP",readable = T)
head(pathway2, 5)

dotplot(pathway2)

write.table(as.data.frame(pathway2),"EnrichedOntology_Pval0.01.txt",sep = "\t",quote = F)
#checking peakAnno and das have same order
tmp = paste(as.data.frame(peakAnno)$seqnames,as.data.frame(peakAnno)$start,sep =":")
tmp = paste(tmp,as.data.frame(peakAnno)$end, sep="-")
if( sum(tmp== rownames(das))== nrow(das)){
  das = cbind(das,as.data.frame(peakAnno))
  write.table(das,"DiffAcc_PVal0.01.txt",sep = "\t",quote = F)
}

Annotation of upregulated peaks only

## Upegulated
das = subset(tt$table,  PValue  < 0.01 & logFC <0)
rownames(das) =paste0("chr",rownames(das))
gr = GRanges(rownames(das))
gr_df = as.data.frame(gr)
peakAnno <- annotatePeak(gr, tssRegion=c(-3000, 3000),TxDb=TxDb.Hsapiens.UCSC.hg38.knownGene, annoDb="org.Hs.eg.db")

>> preparing features information...         2024-01-29 14:12:12 
>> identifying nearest features...       2024-01-29 14:12:12 
>> calculating distance from peak to TSS...  2024-01-29 14:12:12 
>> assigning genomic annotation...       2024-01-29 14:12:12 
>> adding gene annotation...             2024-01-29 14:12:17 

'select()' returned 1:many mapping between keys and columns

>> assigning chromosome lengths          2024-01-29 14:12:17 
>> done...                   2024-01-29 14:12:17 

plotAnnoPie(peakAnno)

gene <- seq2gene(gr, tssRegion = c(-1000, 1000), flankDistance = 3000, TxDb=TxDb.Hsapiens.UCSC.hg38.knownGene)
pathway2 = enrichGO(gene,OrgDb = org.Hs.eg.db,ont="BP",readable = T)
head(pathway2, 5)

dotplot(pathway2)

write.table(as.data.frame(pathway2),"EnrichedOntology_Pval0.01_Up.txt",sep = "\t",quote = F)
#checking peakAnno and das have same order
tmp = paste(as.data.frame(peakAnno)$seqnames,as.data.frame(peakAnno)$start,sep =":")
tmp = paste(tmp,as.data.frame(peakAnno)$end, sep="-")
if( sum(tmp== rownames(das))== nrow(das)){
  das = cbind(das,as.data.frame(peakAnno))
  write.table(das,"DiffAcc_PVal0.01_Up.txt",sep = "\t",quote = F)
}

Annotation of downregulated peaks only

## Downregulated
das = subset(tt$table,  PValue  < 0.01 & logFC >0)
rownames(das) =paste0("chr",rownames(das))
gr = GRanges(rownames(das))
gr_df = as.data.frame(gr)
peakAnno <- annotatePeak(gr, tssRegion=c(-3000, 3000),TxDb=TxDb.Hsapiens.UCSC.hg38.knownGene, annoDb="org.Hs.eg.db")
plotAnnoPie(peakAnno)
gene <- seq2gene(gr, tssRegion = c(-1000, 1000), flankDistance = 3000, TxDb=TxDb.Hsapiens.UCSC.hg38.knownGene)
pathway2 = enrichGO(gene,OrgDb = org.Hs.eg.db,ont="BP",readable = T)
head(pathway2, 5)
dotplot(pathway2)
write.table(as.data.frame(pathway2),"EnrichedOntology_Pval0.01_Down.txt",sep = "\t",quote = F)
#checking peakAnno and das have same order
tmp = paste(as.data.frame(peakAnno)$seqnames,as.data.frame(peakAnno)$start,sep =":")
tmp = paste(tmp,as.data.frame(peakAnno)$end, sep="-")
if( sum(tmp== rownames(das))== nrow(das)){
  das = cbind(das,as.data.frame(peakAnno))
  write.table(das,"DiffAcc_PVal0.01_Down.txt",sep = "\t",quote = F)
}

The files wrote in the text format are available from here

1. Annotation of all peaks

2. Enriched GOs of all peaks

3. Annotation of upregulated peaks

4. Enriched GOs of upregulated peaks

5. Annotation of dowregulated peaks

6. Enriched GOs of downregulated peaks