# RNA-seq Data Processing {#sec-mtd-rna-seq}
The raw RNA-seq datasets (see @tbl-ena below) are available from the ENA website with the following two accessions:
::: {.hover .info #tbl-ena}
| Accession | Description |
|:----------|:------------|
| [PRJEB89714](https://www.ebi.ac.uk/ena/browser/view/PRJEB89714) | ISM2 knock down in EVT(terminal) via small-hairpin RNA |
| [PRJEB110901](https://www.ebi.ac.uk/ena/browser/view/PRJEB110901) | ISM2 knock down in EVT(mid), STB and hTSC via small-hairpin RNA |
RNA-seq datasts deposited in ENA
:::
## Transcript quantification {#sec-mtd-rna-seq-quant}
For an RNA-seq quantification method applied in this study, we used [`Salmon`](https://combine-lab.github.io/salmon/) (v1.10) in mapping-mode to process our RNA-seq datasets.
### Salmon index {#sec-mtd-salmon-index}
```{#lst-salmon-index .bash lst-cap="A pseudo-code to run `salmon index`"}
$HOME/Install/SalmonTools/scripts/generateDecoyTranscriptome.sh \
-j 32 \
-m $HOME/Install/mashmap/mashmap \
-g $HOME/data/genome/Homo_sapiens/Ensembl/GRCh38/Sequence/WholeGenomeFasta/genome.fa \
-a $HOME/data/genome/Homo_sapiens/Ensembl/GRCh38/Annotation/Genes/Homo_sapiens.GRCh38.88.gtf \
-t $HOME/data/genome/Homo_sapiens/Ensembl/GRCh38/Sequence/cDNAFasta/Homo_sapiens.GRCh38.88.cdna.all.ncrna.fa.gz \
-o $HOME/data/Salmon/decoy/Homo_sapiens/Ensembl/GRCh38/Homo_sapiens.GRCh38.88
```
### Salmon quant {#sec-mtd-salmon-quant}
```{#lst-salmon-quant .bash lst-cap="A pseudo-code to run `salmon quant` in mapping mode"}
salmon quant -p 32 \
-i POPS_TR_INDEX \
-l A \
-1 S1_FQ_1.fq \
-2 S1_FQ_2.fq \
--seqBias \
--gcBias \
--posBias \
--discardOrphanQuasi \
--writeUnmappedNames \
--writeMapping \
-o OUTPUT_DIR | samtools view -bS > OUTPUT_DIR/S1_salmon.bam
```
## Differentially expressed gene analysis {#sec-mtd-deg}
### tximeta
The transcript-level read count matrices (e.g. quant.sf files) were imported using [`tximeta`](https://bioconductor.org/packages/release/bioc/html/tximeta.html) (v1.8.5) Bioconductor package and merged at the gene-level.
```{#lst-txi-meta .r lst-cap="R code to construct a transcript-level count matrix via `tximeta`"}
library(magrittr)
library(tximeta)
library(DESeq2)
library(edgeR)
library(BiocParallel)
BiocParallel::register(MulticoreParam(12))
# run info
my.salmon="Salmon"; my.salmon.index="GRCh38.88"
my.slx<-"SLX-23902.Homo_sapiens.v1"
## Tximeta
# Have been done at peta4 - HPCS machine
if(file.exists("RData/ISM2.se.RData")){
load("RData/ISM2.se.RData")
}else{
dt.samples<-data.table(`files`=system(paste0("ls ", "~/results/",my.slx,"/",my.salmon,"/",my.salmon.index,"/*/quant.sf"), intern=T))
dt.samples[,names:=tstrsplit(files,"/",keep=8)]
# Add Condition
dt.samples$Condition<-factor(rep(c("Scr","ISM2"),6),levels=c("Scr","ISM2"))
# Add Pair
dt.samples[,Pair:=paste0("P",ifelse(.I %% 2 ==0,.I/2,(.I+1)/2))]
dt.samples$Pair<-as.factor(dt.samples$Pair)
# check the samples
dt.samples[,-"files"]
# Import quant.sf files
se<-tximeta::tximeta(dt.samples) # tx level
save(se, file="RData/ISM2.se.RData")
}
```
### `DESeq2` {#sec-mtd-deg-deseq2}
We only considered genes found in ≥50% of samples having ≥10 reads, and discarded genes detected as dispersion outliers by DESeq2 (v1.30.0). A total of 17,802 genes and 12 samples were used to find differentially expressed genes by taking the pair information into account in the design matrix of DESeq2. The p-values were calculated from the null hypothesis that the fold changes were less than or equal to 20% (i.e. lfcThreshold=log2(1.2)) in ISM2-shRNA compared to Scr-shRNA (control).
```{r}
#| label: deseq2-setup1
#| eval: false
#| code-summary: "R code to run `DESeq2`"
# filters
minRead=10; minFreq=0.5; minFC=1.2
if(file.exists("RData/ISM2.dds.f2.RData")){
load("RData/ISM2.dds.f2.RData")
}else{
gse<- summarizeToGene(se) # gene-level
# only chr1-22, X, Y, and MT
load("~/data/Annotation/Ensembl/GRCh38.88/dt.ensg.RData")
my.ensg<-dt.ensg[ensembl_gene_id %in% rownames(gse) & chromosome_name %in% c(1:22,"X","Y","MT")]$ensembl_gene_id # %>% length# n=56295
gse<-gse[as.character(my.ensg)]
# set up `dds` at gene-level
my.design <- formula(~ Pair + Condition) # isa 'formula'
dds <- DESeqDataSet(se=gse, design=my.design)
dds <- DESeq(dds, parallel=TRUE) # isa 'DESeqDataSet'
keep <- rowSums(counts(dds) >= minRead) >= ncol(dds)*minFreq # at leat 10 reads for at leat 50 % of the samples
dds.f<- DESeq(dds[keep,], parallel=TRUE) # isa 'DESeqDataSet'
# remove dispOutlier genes
keep2<-!is.na(rowData(dds.f)[,"dispOutlier"]) & !rowData(dds.f)[,"dispOutlier"]
dds.f2<- DESeq(dds.f[keep2], parallel=T)
dim(dds.f2) # 17802 x 12
# just to make sure the base (reference)
dds.f2$Condition<-relevel(dds.f2$Condition, ref="Scr")
save(dds.f2, file="RData/ISM2.dds.f2.RData")
}
my.res<-results(dds.f2,
contrast=c("Condition","ISM2","Scr"), # not necessary for 'ashr'
independentFiltering=FALSE,
lfcThreshold=log2(minFC),
parallel=TRUE)
# apply shink with `res`
my.res3<-lfcShrink(dds.f2,
contrast=c("Condition","ISM2","Scr"), # not necessary for 'ashr'
res=my.res,
lfcThreshold=log2(minFC), # not applicable for 'asher'
type="ashr",
parallel=TRUE)
dt.deseq2<-my.res3 %>% as.data.frame %>% as.data.table(keep.rownames="gene_id")
fwrite(dt.deseq2, file="results/ISM2/ISM2.rna-seq.deseq2.csv")
```
### `edgeR` {#sec-mtd-deg-edger}
For [edgeR](https://bioconductor.org/packages/release/bioc/html/edgeR.html) (v3.32.1) analysis, we used `makeDGEList` function of `tximeta` Bioconductor package to convert the data object of the 17,802 genes across the 12 samples. The gene-level count matrix was normalised by using `calcNormFactors` function of `edgeR` with `TMM` method (trimmed mean of M values) and a quasi-likelihood negative binomial generalised log-linear model (i.e. `glmQLFit`) was applied to account for the batch number, the fetal sex and the gestation information in the design matrix of edgeR. For a statistical test, we used `glmTreat` of edgeR with at least 20% fold-change (i.e. lfc=log2(1.2)).
```{r }
#| label: edger1-setup1
#| eval: false
#| code-summary: "R code to run `edgeR`"
if(!exists("gse")){gse<- summarizeToGene(se)} # gene-level
gse2 <-gse[rownames(dds.f2),]
d2<-tximeta::makeDGEList(gse2)
# TMM normalisation (default). It adds `norm.factors` d2$samples
# NB, we have `offsets`, which take precedence over lib.size and norm.factors
d2<-calcNormFactors(d2,method="TMM")
# design
my.design <- model.matrix(~ 0 +Condition + Pair, data=d2$samples) # isa 'matrix'
my.contrasts <- makeContrasts(`Condition`=ConditionISM2-ConditionScr, levels=my.design)
# dispersion
#dp2 = estimateDisp(d2, design=my.design, verbose=F)
dp2 = estimateGLMCommonDisp(d2, design=my.design, verbose=F)
dp2 = estimateGLMTrendedDisp(dp2, design=my.design, verbose=F)
dp2 = estimateGLMTagwiseDisp(dp2, design=my.design)
# fit
#f = glmFit(dp2, design=my.design) #
f = glmQLFit(dp2, design=my.design) # QL(Quasi-like) pipeline
# get the edgeR results
te <- glmTreat(f, contrast=my.contrasts[,"Condition"], lfc=log2(minFC))
res.edgeR<-topTags(te, n=nrow(te)) # default sort by pvalue
dt.edger<-res.edgeR[["table"]] %>% as.data.table(rownames="gene_id")
# Save DEG from edgeR
fwrite(dt.edger[FDR<0.01], file="results/ISM2/ISM2.rna-seq.DEG.edgeR.csv")
```
## Gene Set Enrichment Analysis
For gene set enrichment analysis (GSEA), we used clusterProfiler (v4.10.1) based on the 17,802 genes ranked in their
descending order using the values from the following formular: sign(log2FC) x -log10(p), where log2FC and p value were from edgeR. The hallmark (collection “H”) and cell type signature (collection “C8”)
gene sets were downloaded using “msigdbr” (v10.0.1) Bioconductor R package and gene sets with a minimum size of 20 genes were considered in GSEA.
```{r}
#| label: gsea-genelist1
#| eval: false
#| code-summary: "R code for Gene Set Enrichment Analysis (GSEA) via `clusterProfiler`"
library(msigdbr)
library(clusterProfiler)
# filter
my.minSize=20
# set ordered gene list
li.geneList<-list()
# gene-set ordered by stat from DESeq2
dt.foo<-merge(dt.deseq2, dt.edger)
dup.genes<-dt.foo[,.N,symbol][N>1]$symbol
dt.deg<-rbind(
dt.foo[symbol %in% dup.genes][order(symbol,-abs(stat))][,.SD[1],symbol],
dt.foo[!symbol %in% dup.genes][order(stat)] %>% setcolorder("symbol")
)
li.geneList[["Pval"]]=dt.deg[order(-sign(logFC)*-log10(PValue))][,.(symbol,foo=sign(logFC)*-log(PValue))]$foo
names(li.geneList[["Pval"]])=dt.deg[order(-sign(logFC)*-log10(PValue))][,.(symbol,foo=sign(logFC)*-log(PValue))]$symbol
# MSig DB
dl.msig<-list()
dl.msig[["Hallmark"]]<-msigdbr(species="human",collection="H") %>% data.table
dl.msig[["Cell Type Signature"]]<-msigdbr(species="human",collection="C8") %>% data.table
li.msig<-
lapply(li.geneList, function(geneList){
lapply(structure(names(dl.msig), names=names(dl.msig)), function(my.msig){
set.seed(20250409)
BPPARAM=BiocParallel::MulticoreParam(workers = 6)
clusterProfiler::GSEA(geneList,
minGSSize = my.minSize, # 10 by default
maxGSSize = 500, # 500 by default
eps = 0, # 1e-50 by default
TERM2GENE = dl.msig[[my.msig]][,.(gs_name,gene_symbol)],
verbose = TRUE,
nPermSimple = 60000, # 1K by default # used by fgsea::fgseaMultilevel
)
})
})
lapply(li.msig, function(li.measure) sapply(li.measure, function(i) nrow(i@result))) # number of significantly enriched terms
dt.MSig<-
lapply(names(li.msig), function(my.measure){
this.set<-sapply(li.msig[[my.measure]], function(i) nrow(i@result)) > 0
lapply(names(li.msig[[my.measure]][this.set]), function(my.msig){
foo<-li.msig[[my.measure]][[my.msig]]
Count=sapply(foo@result$core_enrichment %>% strsplit("/"),length) #%>% unlist #%>% length
GeneRatio=Count / foo@result$setSize
DT<-data.table(Measure=my.measure,Source=my.msig,foo@result, Count, GeneRatio)
DT[,Direction:=ifelse(NES<0,"Down","Up")]
DT[,ID:=gsub("HALLMARK_","",ID)]
}) %>% rbindlist
}) %>% rbindlist
bar<-enrichplot::gseaplot2(li.msig[["stat"]][["Cell Type Signature"]], geneSetID = c("DESCARTES_MAIN_FETAL_EXTRAVILLOUS_TROPHOBLASTS"), pvalue_table = T,base_size=20)
```