test_direct.R

library("IsoformSwitchAnalyzeR")
library(BSgenome)
library(BSgenome.Hsapiens.UCSC.hg38)
library(ggplot2)
Hsapiens


data <- read.csv("/project/owlmayerTemporary/Sid/nanopore-analysis/Direct/Quantification/all_directed_counts.txt")
data <- na.omit(data)
#filter <- read.csv("/project/owlmayerTemporary/Sid/nanopore-analysis/Results_5_1/multiple_exons.txt", header = FALSE)
#data <- data[which(data$transcript %in% filter$V1),]
counts <- data[,c(8,1:6)]
colnames(counts)[1] <- "isoform_id"
counts$isoform_id <- sapply(counts$isoform_id, as.character)

data <- read.csv("/project/owlmayerTemporary/Sid/nanopore-analysis/Direct/Quantification/all_directed_counts_deseq2norm.txt")
data <- na.omit(data)
#filter <- read.csv("/project/owlmayerTemporary/Sid/nanopore-analysis/Results_5_1/multiple_exons.txt", header = FALSE)
#data <- sortdata[which(data$transcript_id %in% filter$V1),]
abundance <- data[,c(9,2:7)]
colnames(abundance)[1] <- "isoform_id"
abundance$isoform_id <- sapply(abundance$isoform_id, as.character)

isoforms <- unique(sort(counts$isoform_id))
write.table(isoforms, file = "/project/owlmayerTemporary/Sid/nanopore-analysis/ReferenceData/isoforms.txt", row.names = FALSE,
            col.names = FALSE, quote = FALSE)

design <- data.frame(sampleID = colnames(abundance)[-1],condition = gsub('_.*', '', colnames(abundance)[-1]))

asswitchlist <- importRdata(isoformExonAnnoation = "/project/owlmayerTemporary/Sid/nanopore-analysis/ReferenceData/gencode.v32.primary_assembly.annotation_tcons.gtf",
                            isoformRepExpression = abundance,
                            isoformCountMatrix = counts, 
                            designMatrix = design,
                            removeNonConvensionalChr = TRUE)


  
  
asswitchlist <- preFilter(switchAnalyzeRlist = asswitchlist,geneExpressionCutoff = 1,isoformExpressionCutoff = 0,removeSingleIsoformGenes = TRUE)
asswitchlist_analyzed <- isoformSwitchTestDEXSeq(switchAnalyzeRlist = asswitchlist, reduceToSwitchingGenes=FALSE)
#asswitchlist_analyzed <- isoformSwitchTestDRIMSeq(switchAnalyzeRlist = asswitchlist, testIntegration = "isoform_only", reduceToSwitchingGenes=TRUE)
asswitchlist_analyzed <- analyzeORF(asswitchlist_analyzed, orfMethod = "longest", genomeObject = Hsapiens)
asswitchlist_analyzed <- extractSequence(asswitchlist_analyzed, pathToOutput = "/project/owlmayerTemporary/Sid/nanopore-analysis/Direct/IsoformAnalyzer/", filterAALength = TRUE, alsoSplitFastaFile = TRUE)


asswitchlist_analyzed <- analyzeAlternativeSplicing(
  switchAnalyzeRlist = asswitchlist_analyzed
)


##################

top <- extractTopSwitches(asswitchlist_analyzed, 
                   filterForConsequences = FALSE, 
                   sortByQvals = TRUE,
                   n = NA)
row.names(top) <- top$Rank

htop <- table(top$gene_switch_q_value)
top <- top[top$gene_switch_q_value %in% names(htop[which(htop == 1)]),] 
sig_genes <- unique(top$gene_name)
write.table(sig_genes[order(sig_genes)], file = "/project/owlmayerTemporary/Sid/nanopore-analysis/Results_5_1/Lists/IsoformSwitchAnalyzeR_List.txt", row.names = FALSE,
            col.names = FALSE, quote = FALSE)


asswitchlist_analyzed

switchPlotTopSwitches(
  switchAnalyzeRlist = asswitchlist_analyzed, 
  n = 10,
  filterForConsequences = TRUE, 
  sortByQvals = TRUE
)


jelena <- read.csv("/project/owlmayerTemporary/Sid/nanopore-analysis/Results_5_1/Lists/JelenaList.txt", header = FALSE)
as.character(jelena$V1) %in% sig_genes 



for (i in 1:5) {
  switchPlot(asswitchlist_analyzed,gene = top$gene_name[i], condition1 = top$condition_1[i], condition2 = top$condition_2[i])
}

mapply(switchPlot,asswitchlist_analyzed,top$gene_name,top$condition_1,top$condition_2)

switchPlot(asswitchlist_analyzed, gene = 'RPS24', condition1 = "day0", condition2 = "day3")
switchPlot(asswitchlist_analyzed, gene = 'RPS24', condition1 = "day0", condition2 = "day5")
switchPlot(asswitchlist_analyzed, gene = 'MYL6', condition1 = "day0", condition2 = "day5")
switchPlot(asswitchlist_analyzed, gene = 'ERBB2', condition1 = "day0", condition2 = "day5")
switchPlot(asswitchlist_analyzed, gene = 'LAGE3', condition1 = "day0", condition2 = "day3")
switchPlot(asswitchlist_analyzed, gene = 'GRWD1', condition1 = "day0", condition2 = "day5")
switchPlot(asswitchlist_analyzed, gene = 'CYTH2', condition1 = "day0", condition2 = "day5")
switchPlot(asswitchlist_analyzed, gene = 'LAD1', condition1 = "day0", condition2 = "day3")
switchPlot(asswitchlist_analyzed, gene = 'EBF2', condition1 = "day0", condition2 = "day3")

switchPlot(asswitchlist_analyzed, gene = 'XLOC_016569', condition1 = "day0", condition2 = "day5")



extractSplicingSummary(
  asswitchlist_analyzed,
  asFractionTotal = FALSE,
  plotGenes=FALSE
)


extractSwitchOverlap(asswitchlist_analyzed)


splicingEnrichment <- extractSplicingEnrichment(
  asswitchlist_analyzed,
  splicingToAnalyze='all',
  returnResult=TRUE,
  returnSummary=TRUE
)

splicingEnrichment <- extractConsequenceEnrichment(
  asswitchlist_analyzed,
  consequencesToAnalyze='all',
  analysisOppositeConsequence=TRUE,
  returnResult=FALSE
)

subset(splicingEnrichment, splicingEnrichment$AStype == 'ATTS')

extractSplicingEnrichmentComparison(
  asswitchlist_analyzed,
  splicingToAnalyze=c('ATTS',"IR"), # Those significant in COAD in the previous analysis
  returnResult=FALSE # Preventing the summary statistics to be returned as a data.frame
)


extractConsequenceSummary(asswitchlist_analyzed, consequencesToAnalyze='all',asFractionTotal = FALSE)


ggplot(data=asswitchlist_analyzed$isoformFeatures, aes(x=dIF, y=-log10(isoform_switch_q_value))) +
  geom_point(
    aes( color=abs(dIF) > 0.3 & isoform_switch_q_value < 0.05 ), # default cutoff
    size=1
  ) +
  geom_hline(yintercept = -log10(0.05), linetype='dashed') + # default cutoff
  geom_vline(xintercept = c(-0.3, 0.3), linetype='dashed') + # default cutoff
  facet_wrap( ~ condition_2) +
  #facet_grid(condition_1 ~ condition_2) + # alternative to facet_wrap if you have overlapping conditions
  scale_color_manual('Signficant\nIsoform Switch', values = c('black','red')) +
  labs(x='dIF', y='-Log10 ( Isoform Switch Q Value )') +
  theme_bw()

######################################

consequencesOfInterest <- c('intron_retention','coding_potential','NMD_status','domains_identified','ORF_seq_similarity')

asswitchlist_test <- analyzeSwitchConsequences(
  asswitchlist_test,
  consequencesToAnalyze = consequencesOfInterest, 
  dIFcutoff = 0.3,
  showProgress=FALSE
)

asswitchlist <- importRdata(isoformExonAnnoation = "/project/owlmayerTemporary/Sid/nanopore-analysis/Results_5_1/GffCompare/nanopore.combined_filt_tcons.gtf",
                            isoformRepExpression = abundance,
                            isoformNtFasta = "/project/owlmayerTemporary/Sid/nanopore-analysis/Results_5_1/Pinfish/corrected_transcriptome_polished_collapsed_tcons.fas",
                            isoformCountMatrix = counts, 
                            designMatrix = design)

exampleSwitchList <- isoformSwitchAnalysisPart1(switchAnalyzeRlist = asswitchlist, genomeObject = Hsapiens, dIFcutoff = 0.3, outputSequences = FALSE, prepareForWebServers = FALSE)

exampleSwitchList <- isoformSwitchAnalysisPart2(switchAnalyzeRlist = exampleSwitchList, dIFcutoff = 0.3, n = 10, removeNoncodinORFs = TRUE,outputPlots = FALSE)

switchPlot(exampleSwitchList,gene = "XLOC_009192", condition1 = "day0", condition2 = "day5")