## ---- include = FALSE----------------------------------------------------
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)

## ---- echo=FALSE, fig.cap="Fig 1. BioTIP workflow with five key analytic steps. RTF: relative transcript fluctuation; PCC: Pearson correlation coefficient; MCI: Module-Criticality Index; Ic: index of critical transition.", fig.align='center', out.width = '60%'----
knitr::include_graphics("Fig1.pdf")

## ------------------------------------------------------------------------
# load package
library(BioTIP)

## ------------------------------------------------------------------------
data(GSE6136_matrix)

dim(GSE6136_matrix)               #[1] 22690rows and 27 columns
row.names(GSE6136_matrix) = GSE6136_matrix$ID_REF
GSE6136 = GSE6136_matrix[,-1]
dim(GSE6136)               #[1] 22690 rows and 26 columns

## ------------------------------------------------------------------------
#requires library(stringr)
library(BioTIP)
data(GSE6136_cli)

#dim(GSE6136_cli) #check dimension
cli = t(GSE6136_cli)

library(stringr)
colnames(cli) = str_split_fixed(cli[1,],'_',2)[,2]
cli = cli[-1,]
cli = data.frame(cli)
cli[,"cell-type:ch1"] = str_split_fixed(cli$characteristics_ch1.1,": ",2)[,2]
cli[,"Ig clonality:ch1"] = str_split_fixed(cli$characteristics_ch1.3,": ",2)[,2]

colnames(cli)[colnames(cli) == "cell-type:ch1"] = "group"
cli$Row.names = cli[,1]
head(cli[,1:3])

## ------------------------------------------------------------------------
dat <- GSE6136
df <- log2(dat+1)
head(df)

## ---- warning=FALSE------------------------------------------------------
tmp <- names(table(cli$group))
samplesL <- split(cli[,1],f = cli$group)
#head(samplesL)
test <- sd_selection(df, samplesL,0.01)
head(test[["activated"]])

## ----echo=TRUE, warning=FALSE--------------------------------------------
library(BioTIP)
#library(igraph)
library(cluster)
igraphL <- getNetwork(test, fdr = 1)

cluster <- getCluster_methods(igraphL)

## ----echo=TRUE, warning=FALSE--------------------------------------------
names(cluster)
head(cluster[[1]])

## ----echo=TRUE, warning=FALSE--------------------------------------------
membersL_noweight <- getMCI(cluster,test,adjust.size = FALSE)
plotBar_MCI(membersL_noweight,ylim = c(0,6))

## ----echo=TRUE, warning=FALSE--------------------------------------------
maxMCIms <- getMaxMCImember(membersL_noweight[[1]],membersL_noweight[[2]],min =10)
names(maxMCIms)
names(maxMCIms[[1]])
names(maxMCIms[[2]])

## ----echo=TRUE, warning=FALSE--------------------------------------------
head(maxMCIms[['idx']])
head(maxMCIms[['members']][['lymphoma_aggressive']])

## ------------------------------------------------------------------------
biomodules = getMaxStats(membersL_noweight[['members']],maxMCIms[[1]])
maxMCI = getMaxStats(membersL_noweight[['MCI']],maxMCIms[[1]])
head(maxMCI)

## ------------------------------------------------------------------------
maxSD = getMaxStats(membersL_noweight[['sd']],maxMCIms[[1]])
head(maxSD)

## ------------------------------------------------------------------------
CTS = getCTS(maxMCI, maxMCIms[[2]])

## ----echo=TRUE, warning=FALSE--------------------------------------------
par(mar = c(10,5,0,2))
plotMaxMCI(maxMCIms,membersL_noweight[[2]],states = names(samplesL),las = 2)

## ----echo=TRUE, warning=FALSE--------------------------------------------
simuMCI <- simulationMCI(3,samplesL,df)
plot_MCI_Simulation(maxMCI,simuMCI)

## ------------------------------------------------------------------------
IC <- getIc(df,samplesL,CTS)
par(mar = c(10,5,0,2))
plotIc(IC,las = 2)

## ----warning=FALSE-------------------------------------------------------
simuIC <- simulation_Ic(length(CTS),samplesL,df)
par(mar = c(10,5,0,2))
plot_Ic_Simulation(IC,simuIC,las = 2)

## ------------------------------------------------------------------------
simulated_Ic = plot_simulation_sample(df,length(samplesL[['lymphoma_aggressive']]),IC[['lymphoma_aggressive']],CTS)

## ---- echo=FALSE, fig.align='center', out.width = '65%'------------------
knitr::include_graphics("Fig2.jpg")

## ------------------------------------------------------------------------
library(BioTIP)

data(gencode)
head(gencode)

## ----"python code", eval = FALSE-----------------------------------------
#  gtf = ("Your/PATH/TO/THE/FILE")
#  outF = open("gtf_summary_transbiotype.txt","w")
#  
#  def getquote(str,f,target):
#      targetLen = len(target)+2
#      strInd = str.find(target)
#      st = strInd + len(target)+2
#      ed = st + str[st:].find("";")
#      #print(st,ed)
#      f.write("\t"+str[st:ed]) if strInd!= -1 else f.write("\t"+"NA.")
#  
#  with open(gtf, "r") as f:
#       for line in f:
#          if line[0] != "#":
#              chromosome = line.split("\t")[0]
#              st = line.split("\t")[3]
#              ed = line.split("\t")[4]
#              strand = line.split("\t")[6]
#              type = line.split("\t")[2]
#              outF.write(chromosome+"\t"+st+"\t"+ed+"\t"+strand+"\t"+type)
#              c = "transcript_id"
#              g = "gene_name"
#              t = "transcript_type"
#              getquote(line,outF,c)
#              getquote(line,outF,g)
#              getquote(line,outF,t)
#              outF.write("\n")
#  outF.close()

## ------------------------------------------------------------------------
library(BioTIP)
library(GenomicRanges)
data(gencode)
head(gencode)

## ------------------------------------------------------------------------
query <- GRanges(c("chr1:2-10:+","chr1:6-10:-"), Row.names = c("trans1","trans2"), score = c(1,2))
head(query)


## ------------------------------------------------------------------------
library(BioTIP)
gr <- GRanges(c("chr1:1-5:+","chr1:2-3:+"),biotype = c("lincRNA","CPC"))
head(gr)

## ------------------------------------------------------------------------

intron <- GRanges("chr1:6-8:+")
head(intron)

## ------------------------------------------------------------------------

intron_trncp <- getBiotypes(query, gr, intron)
intron_trncp

## ------------------------------------------------------------------------
library(BioTIP)
data("intron")
data("ILEF")
data("gencode")

gencode_gr = GRanges(gencode)
ILEF_gr = GRanges(ILEF)
cod_gr = GRanges(cod)
intron_gr = GRanges(intron)

non_coding <- getBiotypes(ILEF_gr, gencode_gr, intron_gr)
dim(non_coding)
head(non_coding[,1:3])

## ------------------------------------------------------------------------
coding <- getBiotypes(ILEF_gr, gencode_gr)
dim(coding)
head(coding[,1:3])

## ------------------------------------------------------------------------
library(BioTIP)

data(ILEF)
data(cod)
ILEF_gr = GRanges(ILEF)
cod_gr = GRanges(cod)

rdthrough <- getReadthrough(ILEF_gr, cod_gr)
head(rdthrough)

## ----SessionInfo---------------------------------------------------------
sessionInfo()