RSVinfected Host cells

Count reads map to hg38 genome

library(m6Amonster)
samplenames <-c("Ctl","RSVinfect1","RSVinfect2")
gtf <- "~/Database/genome/hg38/hg38_UCSC.gtf"
RSV.A549 <- countReads(samplenames = paste(samplenames,"align_hg38",sep = "." ),
                          gtf = gtf,
                          bamFolder = "/home/zijiezhang/RSV/201803/bam_files",
                          outputDir = "/home/zijiezhang/RSV/201803",
                          modification = "m6A",
                          threads = 40
                          )
save(RSV.A549,file = "~/RSV/201803/RSV.A549_infected.RData")

Call peaks

library(m6Amonster)
load("~/RSV/201803/RSV.A549_infected.RData")
#RSV.A549 <- callPeakFisher(RSV.A549,threads = 20)

Report peaks for Ctl and Infected samples

CtlPeak <- reportConsistentPeak(RSV.A549,samplenames = "Ctl.align_hg38",threads = 15)

Reporting peak concsistent in all samples for
 Ctl.align_hg38 
Hyper-thread registered: TRUE 
Using 15 thread(s) to report merged report...
Time used to report peaks: 2.7153252919515 mins... 

InfectedPeak <- reportConsistentPeak(RSV.A549,samplenames = c("RSVinfect1.align_hg38","RSVinfect2.align_hg38"), threads = 15 )

Reporting peak concsistent in all samples for
 RSVinfect1.align_hg38 RSVinfect2.align_hg38 
Hyper-thread registered: TRUE 
Using 15 thread(s) to report merged report...
Time used to report peaks: 1.45808432896932 mins... 

Plot distribution of peaks on the meta-gene

plotMetaGeneMulti(list("Ctl"=CtlPeak, "RSV infected" = InfectedPeak), "~/Database/genome/hg38/hg38_UCSC.gtf")

[1] "Converting BED12 to GRangesList"
[1] "It may take a few minutes"
[1] "Converting BED12 to GRangesList"
[1] "It may take a few minutes"
[1] "total 58259 transcripts extracted ..."
[1] "total 42543 transcripts left after ambiguity filter ..."
[1] "total 21293 mRNAs left after component length filter ..."
[1] "total 7986 ncRNAs left after ncRNA length filter ..."
[1] "Building Guitar Coordinates. It may take a few minutes ..."
[1] "Guitar Coordinates Built ..."
[1] "Using provided Guitar Coordinates"
[1] "resolving ambiguious features ..."
[1] "no figure saved ..."

Please note: This plot only shows m6A peak distribution of both Ctl and infected samples are consistent with conanical pattern, suggesting the m6A-seq is of okay quality. It doesn’t tell any difference between the two group! The peak show up in one group but not the other group are just some boundary cases where p-value are close to the threshold for peak-calling. Besides, Ctl has no replicate while Infected peaks are consistent peaks of two replicates(more stringent); resulting in much fewer peaks in Infected group.

Show Ctl peaks distribution on Pie chart

peakDistribution(peak = CtlPeak, gtf = "~/Database/genome/hg38/hg38_UCSC.gtf" )

Show RSV infected peaks distribution on Pie chart

peakDistribution(peak = InfectedPeak, gtf = "~/Database/genome/hg38/hg38_UCSC.gtf" )

Differential methylation test with QNB

#RSV.A549 <- JointPeakCount(RSV.A549,joint_threshold = 2)
#RSV.A549 <- normalizePeak(RSV.A549,c("Ctl","Infected","Infected"))
#RSV.A549 <- QNBtest(RSV.A549)
hist(qvalue::qvalue(RSV.A549$all.est.peak$p_value))

Display differential methylated peaks

jointPeak <- reportJointPeak(RSV.A549,threads = 10)

Reporting joint peak with PoissonGamma test statistics...
Reporting  16156  peaks...
Hyper-thread registered: TRUE 
Using 10 thread(s) to report merged report...
Time used to report peaks: 4.10002266565959 mins... 

#write.table(jointPeak[which(jointPeak$fdr<0.1 & abs(jointPeak$logFC)>0.5 ),], file = "~/RSV/201803/Host_diff_meth.xls",sep = "\t", quote = F,row.names = F)
knitr::kable(jointPeak[which(jointPeak$fdr<0.1 & abs(jointPeak$logFC)>0.5 ),], format = "pandoc", row.names = F)

** Please note: ** there are very few differentially methylated sites on the host cells in this data set. Those peaks are ones that are different between Ctl and Infected group.

Differential expression analysis

The “emerged/lost peak figure (analysis)” before was an incorrect way to analyzing the data. Those emerged/lost peads mostly reflected differentially expressed genes.

deseq2_res <- RADAR::DESeq2(RSV.A549$geneSum, matrix(c("Ctl","Infected","Infected"),ncol = 1) )

[1] "without pc ..."
Using input matrix without additional normalization...

#write.table(deseq2_res[which(deseq2_res$p.adjust<0.05),], file = "~/RSV/201803/Host_diff_expression.xls",sep = "\t", quote = F,row.names = F)
knitr::kable(deseq2_res[which(deseq2_res$p.adjust<0.1),], format = "pandoc", row.names = F)