1

arXiv:1403.5029v3 [] 15 Sep 2015

Network-based Isoform Quantification with RNA-Seq Data for
Cancer Transcriptome Analysis
Wei Zhang1 , Jae-Woong Chang2 , Lilong Lin3 , Kay Minn4 , Baolin Wu5 , Jeremy Chien4 , Jeongsik Yong2 , Hui
Zheng3 , and Rui Kuang1,∗
1 Department of Computer Science and Engineering, University of Minnesota Twin Cities,
Minneapolis, Minnesota, United States of America
2 Department of Biochemistry, Molecular Biology and Biophysics, University of Minnesota Twin
Cities, Minneapolis, Minnesota, United States of America
3 Guangzhou Institutes of Biomedicine and Health, Chinese Academy of Sciences, Guangzhou,
Peoples Republic of China
4 Department of Cancer Biology, University of Kansas Medical Center, Kansas City, Kansas,
United States of America
5 Division of Biostatistics, School of Public Health, University of Minnesota Twin Cities,
Minneapolis, Minnesota, United States of America
∗ E-mail: kuang@cs.umn.edu

Abstract
High-throughput mRNA sequencing (RNA-Seq) is widely used for transcript quantification of gene isoforms.
Since RNA-Seq data alone is often not sufficient to accurately identify the read origins from the isoforms for
quantification, we propose to explore protein domain-domain interactions as prior knowledge for integrative analysis with RNA-Seq data. We introduce a Network-based method for RNA-Seq-based Transcript Quantification
(Net-RSTQ) to integrate protein domain-domain interaction network with short read alignments for transcript
abundance estimation. Based on our observation that the abundances of the neighboring isoforms by domaindomain interactions in the network are positively correlated, Net-RSTQ models the expression of the neighboring
transcripts as Dirichlet priors on the likelihood of the observed read alignments against the transcripts in one
gene. The transcript abundances of all the genes are then jointly estimated with alternating optimization of
multiple EM problems. In simulation Net-RSTQ effectively improved isoform transcript quantifications when
isoform co-expressions correlate with their interactions. qRT-PCR results on 25 multi-isoform genes in a stem
cell line, an ovarian cancer cell line, and a breast cancer cell line also showed that Net-RSTQ estimated more
consistent isoform proportions with RNA-Seq data. In the experiments on the RNA-Seq data in The Cancer
Genome Atlas (TCGA), the transcript abundances estimated by Net-RSTQ are more informative for patient
sample classification of ovarian cancer, breast cancer and lung cancer. All experimental results collectively
support that Net-RSTQ is a promising approach for isoform quantification. Net-RSTQ toolbox is available at
http://compbio.cs.umn.edu/Net-RSTQ/.

Author Summary
New sequencing technologies for transcriptome-wide profiling of RNAs have greatly promoted the interest in
isoform-based functional characterizations of a cellular system. Elucidation of gene expressions at the isoform
resolution could lead to new molecular mechanisms such as gene-regulations and alternative splicings, and potentially better molecular signals for phenotype predictions. However, it could be overly optimistic to derive the
proportion of the isoforms of a gene solely based on short read alignments. Inherently, systematical sampling
biases from RNA library preparation and ambiguity of read origins in overlapping isoforms pose a problem in
reliability. The work in this paper exams the possibility of using protein domain-domain interactions as prior
knowledge in isoform transcript quantification. We first made the observation that protein domain-domain interactions positively correlate with isoform co-expressions in TCGA data and then designed a probabilistic EM
approach to integrate domain-domain interactions with short read alignments for estimation of isoform propor-

2

tions. Validated by qRT-PCR experiments on three cell lines, simulations and classifications of TCGA patient
samples in several cancer types, Net-RSTQ is proven a useful tool for isoform-based analysis in functional genomes
and systems biology.

Introduction
Application of next generation sequencing technologies to mRNA sequencing (RNA-Seq) is a widely used approach
in transcriptome study [1, 2, 3]. Compared with microarray technologies, RNA-Seq provides information for
expression analysis at transcript level and avoids the limitations of cross-hybridization and restricted range of
the measured expression levels. Thus, RNA-Seq is particularly useful for quantification of isoform transcript
expressions and identification of novel isoforms. Accurate RNA-Seq-based transcript quantification is a crucial
step in other downstream transcriptome analyses such as isoform function prediction in the pioneer work in
[4], and differential gene expression analysis [5] or transcript expression analysis [6]. Detecting biomarkers from
transcript quantifications by RNA-Seq is also a frequent common practice in biomedical research. However,
transcript quantification is challenging since a variety of systematical sampling biases have been observed in
RNA-Seq data as a result of library preparation protocols [7, 8, 9, 10]. Moreover, in the aligned RNA-Seq short
reads, most reads mapped to a gene are potentially originated by more than one transcript. The ambiguous
mapping could result in hardly identifiable patterns of transcript variants [10, 11].
A useful prior knowledge that has been largely ignored in RNA-Seq transcriptome quantification is the relation
among the isoform transcripts by the interactions between their protein products. The protein products of
different isoforms coded by the same gene may contain different domains interacting with the protein products
of the transcripts in other genes. Previous studies suggested that alternative splicing events tend to insert or
delete complete protein domains/functional motifs [12] to mediate key linkages in protein interaction networks by
removal of protein domain-domain interactions [13]. The work in [4, 14] also suggested unique patterns in isoform
co-expressions. Thus, the abundance of an isoform transcript in a gene can significantly impact the quantification
of the transcripts in other genes when their protein products interact with each other to accomplish a common
function as illustrated by a real subnetwork in Fig 1, which is constructed based on domain-domain interaction
databases [15, 16] and Pfam [17]. Motivated by our observation that the protein products of highly co-expressed
transcripts are more likely to interact with each other by protein domain-domain binding in four TCGA RNA-Seq
datasets (see the section Results), we constructed two human transcript interaction networks of different sizes
based on protein domain-domain interactions to improve transcript quantification. Based on the constructed
transcript network, we propose a network-based transcript quantification model called Net-RSTQ to explore
domain-domain interaction information for estimating transcript abundance. In the Net-RSTQ model, Dirichlet
prior representing prior information in the transcript interaction network is introduced into the likelihood function
of observing the short read alignments. The new likelihood function of Net-RSTQ can be alternating-optimized
over each gene with expectation maximization (EM). It is important to note that the Dirichlet prior from the
neighboring isoforms play two possible roles. On one hand, for the isoforms in the same gene but with different
interacting partners, the different prior information will help differentiate their expressions to reflect their different
functional roles. On the other hand, for the isoforms in the same gene with the same interacting partners, the
uniform prior assumes no difference in their functional roles and thus, promotes a smoother expression patterns
across the isoforms. In both cases, the Dirichlet prior captures the functional variations/similarities across the
isoforms in each gene as prior information for estimation of their abundance.

3

Figure 1. An isoform transcript network based on protein domain-domain interactions. (A) The
subnetwork shows the domain-domain interactions among transcripts from four human genes, CD79B, CD79A,
LCK and SYK. In the network, the nodes represent isoform transcripts, which are further grouped and
annotated by their gene name; and the edges represent domain-domain interactions between two transcripts.
Each edge is also annotated by the interacting domains in the two transcripts. (B) RefSeq transcript
annotations of CD79A and CD79B are shown with Pfam domain marked in color. The Pfam domains were
detected with Pfam-Scan software. Note that no interaction is included between transcripts NM 001039933 and
NM 000626 of gene CD79B without assuming self-interactions for modeling simplicity. For better visualization,
only the interactions coincide with PPI are shown in the figure.
The paper is organized as following. In the section Materials and Methods, we describe the procedure to
construct protein domain-domain interaction networks, the mathematic description of the probabilistic model and
the Net-RSTQ algorithm, qRT-PCR experiment design, and RNA-Seq data preparation. In the section Results,
we first demonstrate the correlation between protein domain-domain interactions and isoform transcript coexpressions across samples in four cancer RNA-Seq datasets from The Cancer Genome Atlas (TCGA) to justify
using domain-domain interactions as prior knowledge. We then compared the predicted isoform proportions
with qRT-PCR experiments on 25 multi-isoform genes in three cell lines, H9 stem cell line, OVCAR8 ovarian
cancer cell line and MCF7 breast cancer cell line. Net-RSTQ was also applied to four cancer RNA-Seq datasets
to quantify isoform expressions to classify patient samples by the survival or relapse outcomes. In addition,
simulations were also performed to measure the statistical robustness of Net-RSTQ over randomized networks.

Materials and Methods
In this section, we first describe the construction of the transcript interaction network and review the base
probabilistic model for transcript quantification with RNA-Seq data. We then introduce the network-based
transcript quantification model (Net-RSTQ) by applying the protein domain-domain interaction information as
prior knowledge. The notations used in the equations are summarized in Table 1. At last, qRT-PCR experiment
design and RNA-Seq data preparation are explained.

4
Notation
N
T
lik
r
pik
Pi
P
ρik
π
φik
α
qijk
S

Description
total # of genes
set of transcripts; Tik is the kth transcript of the ith gene; Ti denotes the transcripts of the ith gene
length of transcript Tik
set of reads; rij is the j th read aligned to the ith gene; ri is the read set aligned to the ith gene
the probability of a read generated by transcript Tik in the ith gene
the probability of a read generated by transcript Ti in the ith gene, specifically, [pi1 , ..., pi,|Ti | ]
concatenate of all Pi , specifically, [P1 , P2 , ..., PN ]
relative abundance of the transcript Tik in the ith gene
transcript expression; πik is the expression of the kth transcript of the ith gene
average expressions (normalized) of transcript Tik ’s neighbors in the transcript network
parameters of Dirichlet distribution; αik = λφik + 1 is the parameter of the Dirichlet distribution of pik
read sampling probability, qijk = lik −l1 r +1 if read rij is aligned to transcript Tik , otherwise qijk = 0
binary matrix for transcript interaction network

Table 1. Notations

Transcript network construction
Two binary transcript networks were constructed by measuring the protein domain-domain interactions (DDI)
between the domains in each pair of transcripts in four steps. First, the translated transcript sequences of all
human genes were obtained from RefSeq [18]. Second, Pfam-Scan was used to search Pfam databases for the
matched Pfam domains on each transcript with 1e-5 e-value cutoff [17]. Note that only high quality, manually
curated Pfam-A entries in the database were used in the search. Third, domain-domain interactions were obtained
from several domain-domain interaction databases, and if any domain-domain interaction exists between a pair of
transcripts, the two transcripts are connected in the transcript network. Specifically, 6634 interactions between
4346 Pfam domain families from two 3D structure-based DDI datasets (iPfam [15] and 3did [16]) inferred from
the protein structures in Protein Data Bank (PDB) [19] were used in the experiments. Besides these highly
confident structure-based DDIs, transcript interactions constructed from 2989 predicted high-confidence DDIs
and 2537 predicted medium-confidence DDIs in DOMINE [20] were also included if the transcript interaction
agrees with protein-protein interactions (PPI) in HPRD [21].
In the experiments, we focused on the transcripts from two cancer gene lists from the literature for better
reliability in annotations. The first smaller transcript network consists of 11736 interactions constructed from the
3D structure-based DDIs and 421 interactions constructed from the predicted DDIs among the 898 transcripts
in 397 genes from the first gene list [22]. The second larger transcript network contains 711,516 interactions
constructed from the 3D structure-based DDIs among 5599 transcripts in 2551 genes in a larger gene list [23].
Since inclusion of the predicted DDIs results in a much higher density in the large network, the large network
does not include predicted DDIs to prevent too many potential false positive interactions. The characteristics of
the two transcript networks are summarized in Table 2. The density of the two networks are 3.02% and 4.54%
respectively, which are in similar scale with the PPI network. Both networks show high clustering coefficients,
suggesting modularity of subnetworks. Note that self-interactions (interactions between transcript(s) in the
same gene) are not considered since Net-RSTQ only utilizes positive correlation between the expressions of
neighboring transcripts in different genes. For simplicity, Net-RSTQ assumes that self-interactions will not
change the transcript quantification of an individual gene in the model.
Small Network
Large Network

# of Gene
397
2551

# of Transcripts
898
5599

# of Interactions
12157
711516

Density
3.02%
4.54%

Diameter
9
9

Avg. # of Neighbors
27.08
254.16

Avg. Cluster Coefficients
0.3578
0.5255

Table 2. Network characteristics.
In Fig 1(A) a subnetwork of the transcripts in gene CD79A and CD79B with their direct neighbors in the small

5

transcript network is shown. The RefSeq transcript annotations of CD79A and CD79B are shown in Fig 1(B). In
CD79A transcript NM 001783 contains an extra domain pfam07686 while transcript NM 021601 only contains a
shorter hit pfam02189. Note pfam02189 also has the same hit in NM 001783 with an e-value larger than 1e-5. In
CD79B transcripts NM 001039933 and NM 000626 contain a domain pfam07686, which is removed in alternative
splicing of NM 021602. In the transcript subnetwork shown in Fig 1(A), the transcripts in CD79A or CD79B
have different interaction partners in the network. In the transcripts in CD79A, the expression of NM 021601 will
correlate with the transcripts in LCK and SYK, and NM 001783 will correlate with two transcripts in CD79B.
The isoform transcripts in LCK and SYK show no different DDIs suggesting there is no functional variation by
protein bindings and more similar expression patterns are potentially expected as prior knowledge.

Base model for transcript quantification
We first consider the method proposed in [24, 25] as the base model for quantification of the transcripts in a
single gene. Let Ti denote the set of the transcripts in the ith gene and Tik be the kth transcript in Ti . The
probability of a read being generated by the transcripts in Ti is modeled by a categorical distribution specified
P|Ti |
by parameters pik , where k=1
pik = 1 and 0 ≤ pik ≤ 1. For the set of the reads ri aligned to gene i, we
consider the likelihood of that each of the |ri | short reads is sampled from one of the transcripts to which the
read aligns. Specifically, for each read rij aligned to transcript Tik , the probability of obtaining rij by sampling
from Tik , namely P r(rij |Tik ) is qijk = lik −l1r +1 [26, 8, 27], where lr is the length of the read. Assuming each
read is independently sampled from one transcript, the uncommitted likelihood function [24] to estimate the
parameters Pi from the observed read alignments against gene i is
L(Pi ; ri ) = P r(ri |Pi )

=

|ri |
Y

P r(rij |Pi ) =

j=1

|ri | |Ti |
Y
X

P r(Tik |Pi )P r(rij |Tik ) =

j=1 k=1

|ri | |Ti |
Y
X

pik qijk .

(1)

j=1 k=1

This likelihood function is concave but it may contain plateau in the likelihood surface. Therefore, Expectation
Maximization (EM) is then applied to obtain the optimal Pi . In the EM algorithm, the expectation of read
assignments to transcripts were estimated in the E-step and the likelihood function with the expected assignments
can be maximized in the M-step to estimate Pi . The relative abundance of the transcript Tik in gene i, ρik , can
be derived from
pik

ρik = P|Tlik|
i

pik
k=1 lik

,

(2)

and the transcript expressions in gene i, πik , can be calculated by
πik =

|ri |pik
.
lik

(3)

The base model is applied independently to each individual gene and no relation among the transcripts is
considered.

Network-based transcript quantification model
In the Net-RSTQ model, the transcript interaction network S based on protein domain-domain interactions is
introduced to calculate a prior distribution for estimating P jointly across all the genes and all the transcripts.
The model assumes that the prior distribution of Pi is a Dirichlet distribution specified by parameters αi and
each αik is proportional to the read count by average expression of the transcript Tik ’s neighbors in the transcript
network S. The prior read count φik is defined as follows,

6

S∗,(i,k)
0
φik = lik (π P
),
(S∗,(i,k) )

(4)

where S∗,(i,k) is a binary vector represents the neighborhood of transcript Tik in transcript network S and
P
(S∗,(i,k) ) is the size of the neighborhood. The calculation of each φik is illustrated in Fig 2. The Dirichlet
parameter αi is defined as a function of φik as
αik = λφik + 1,

(5)

where λ > 0 is a tuning parameter balancing the belief between the prior-read count and the aligned-read count.

Figure 2. Transcript interaction neighborhood. In this toy example, transcript Tik has four neighbor
transcripts {Tg1 a , Tg2 b , Tg3 c , Tg4 d }, which are transcripts from g1 , g2 , g3 and g4 , respectively. The neighborhood
expression φik of Tik is then calculated as the average of its neighbor transcripts’ expressions and further
normalized by transcript
P length, represented as the vector product between π and S∗,(i,k) normalized by the
number of neighbors
S∗,(i,k) and the transcript length lik in the figure.
To obtain the optimal P jointly for all genes, we introduce a pseudo-likelihood model to estimate P iteratively
in each iteration. Assuming uniform P r(ri ), the pseudo-likelihood function is defined as,
L(P , α; r) =

N
Y

L(Pi , αi ; ri ) =

N
Y
P r(Pi |αi )P r(ri |Pi )

P r(ri )

i=1

i=1

∝

N
Y

P r(Pi |αi )P r(ri |Pi ).

(6)

i=1

Note that the pseudo-likelihood model relies on the independence assumption among the likelihood functions of
each individual gene when the α parameters of the Dirichlet priors are pre-computed. Thus, the model simply
takes the product of the likelihood function from each gene. Each prior distribution P r(Pi |αi ) follows the
Dirichlet distribution,
|Ti |

P r(Pi |αi ) = C(αi )

Y
k=1

P
Γ(
αik )
pik αik −1 , where C(αi ) = Q k
.
k Γ(αik )

(7)

Integrating equations (1) and (7), the pseudo-likelihood function in equation (6) can be rewritten with Dirichlet prior as

L(P ; r)

=

N
Y



|Ti |



C(αi )

Y

pik αik −1  

i=1

=

N
Y
i=1

|ri | |Ti |
Y
X


pik qijk 

j=1 k=1

k=1



|Ti |



C(λφi + 1)

Y

pik λφik  

k=1

|ri | |Ti |
Y
X

j=1 k=1


pik qijk  .

(8)

7

In the pseudo-likelihood function in equation (8), the only hyper-parameter λ balances the proportion between
the Dirichlet priors and the observed read counts of each transcript. The larger the λ, the more belief put on the
priors.

The Net-RSTQ algorithm
The Net-RSTQ algorithm optimizes equation (8) by dividing the optimization into sub-optimization problems of
sequentially estimating each Pi . Specifically, we fix all Pc , c 6= i, and thus φi when estimating Pi with EM in
each iteration and repeat the process multiple rounds throughout all the genes. In each step, the neighborhood
expression φ is recomputed with new Pi for computing the quantification of the next gene. For each suboptimization problem, we estimate Pi with a fixed φ, the part of the likelihood function in equation (8) involved
with the current variables Pi is




|Tg |
|ri | |Ti |
|Ti |
Y
X
Y
Y
Y
(9)
pik qijk  ,
pik λφik  
L̄(Pi ; ri ) = 
C(λφg + 1)
pgk λφgk  C(λφi + 1)
g∈nb(i)

k=1

k=1

j=1 k=1

where nb(i) is the set of the genes containing transcripts that are neighbors of the transcripts in gene i in the
transcript network. Equation (9) consists of three terms separated by the braces. The second and the third
terms are the Dirichlet prior and the likelihood of the observed counts in the data for gene i. The first term
is the Dirichlet priors of the neighbor transcripts of each Tik . These prior probabilities are involved since φg
are functions of the current variable Pi (equations (3)-(5)). Equation (9) cannot be easily solved with standard
techniques. We adopt a heuristic approach to only take steps that will increase the whole pseudo-likelihood
function in equation (8). The Net-RSTQ algorithm is outlined below
Algorithm 1 Net-RSTQ
1: Initialization: random initialization or base EM (equation (1)) estimation of P (0)
2: for round t = 1, . . . do
3:
P (t) = P (t−1)
4:
for gene i = 1, . . . , N do
5:
compute φi based on P (t) with equations (3) and (4)
6:
estimate Pi with EM algorithm (see next section)
(t)

7:
if L̄(Pi ) > L̄(Pi ) then
8:
Pi (t) = Pi
9:
end if
10:
end for
11:
if max(abs(P (t) − P (t−1) )) <1e-6 then
12:
break
13:
end if
14: end for
15: return P

In the algorithm, the outer for-loop between line 2-14 performs multiple passes of updating P . The inner
for-loop between line 4-10 scans through each gene to update each Pi . Line 7 checks the the difference in the
likelihood L̄ of gene i before and after the estimated Pi is applied. The newly estimated Pi is kept in line
8 only if the likelihood L̄ in equation (9) is higher. The convergence of P is checked at line 11. In each suboptimization problem, EM algorithm (described in the next section) is applied to estimate Pi . After convergence,
the transcripts expression π can be learned by equation (3) with the optimal P .

8

Estimating Pi given φi
In line 6 of Algorithm 1, we maximize the likelihood function of the sub-optimization problem in equation (9) to
learn Pi as



|Ti |
|ri | |Ti |
Y
Y
X
L(Pi ; ri ) = C(λφi + 1)
pik λφik  
pik qijk  .
(10)
j=1 k=1

k=1

Note that equation (10) is the part of equation (9) without the Dirichlet priors of the neighboring genes. In line
7 of Algorithm 1, the ignored Dirichlet priors are combined with the likelihood in equation (10), when L̄(Pi ) is
computed, to evaluate the whole likelihood in equation (9). The likelihood function in equation (10) is defined on
a categorical variable with Dirichlet prior, which can be solved with EM algorithm. Following EM formulation in
[26], the expectation aijk , a soft assignment of read j to transcript k in gene i, is first estimated in the expectation
step and Pi is then learned in the maximization step. When φi is given, by taking log of equation (10) we can
write the EM steps to find Pi below.
E step:
Letting Match signify a matching between reads and transcripts, and Match(j) be the transcript from which
read j originates, we get:
log[L(Pi ; ri , M atch)]

=

logC(λφi + 1) +

|Ti |
X

λφik log(pik ) +

|ri |
X

log(piM atch(j) qijM atch(j) ),

(11)

j=1

k=1

which leads to
(it)

Q(Pi |Pi

)

= EM atch|ri ,P (it) [log(L(Pi ; ri ))]
i

=

logC(λφi + 1) +

|Ti |
X

(it)

λφik log(pik ) +

k=1

=

logC(λφi + 1) +

|Ti |
X

(it)
λφik log(pik )

+

|ri | |Ti |
X
X

(it)

p qijk
(it)
(logpik + logqijk ) ∗ P|T ik| (it)
i
j=1 k=1
k=1 pik qijk

|ri | |Ti |
X
X

(it)
aijk log(pik )

+

j=1 k=1

k=1

|ri | |Ti |
X
X

aijk log(qijk )

(12)

j=1 k=1

where it is the itth iteration in EM and
(it)

p qijk
aijk = P|T ik| (it)
.
i
k=1 pik qijk
M step:
Given that qijk and φi are known, the above reduces to maximizing


|Ti |
|ri | |Ti |
X
X
X
(it+1)
Pi
= arg max 
λφik log(pik ) +
aijk log(pik ) .
Pi

k=1

(13)

(14)

j=1 k=1

Using Lagrange multipliers and differentiating, equation (14) is maximized when
P|ri |
aijk
λφik + j=1
(it+1)
.
pik
= P|T |
P|ri |
i
j=1 aijk )
k=1 (λφik +

(15)

After EM algorithm converges, we update P with the newly estimated Pi only if the update leads to increase of
equation (9). It can be seen from equation (15) that the role of λ is a parameter controlling the balance between

9

the prior-read count and the aligned-read count. To see that, recall φik is the prior-read count of transcript Tik
P|ri |
by the average expression of its neighbors (equation (4)) and j=1
aijk is the expected aligned-read count of
transcript Tik . λ directly balances the contributions from the two terms. Therefore, a reasonable choice of λ
should apply to RNA-Seq data with similar level of noise or bias in general.

qRT-PCR experiment design
Three qRT-PCR experiments are designed to measure the isoform proportions of 25 multi-isoform genes in
three cell lines, H9 stem cell line, OVCAR8 ovarian cancer cell line and MCF7 breast cancer cell line. The cell
lines were selected based on the available of both RNA-Seq data and cell culture in our labs. The qRT-PCR
experiments focused on the gene with most different quantification results reported by Net-RSTQ and other
compared methods. Due to the limitations in time and cost of running qRT-PCR experiments, only the 25 genes
in the three cell lines were tested with all the results reported in the experiments. Quantitation of the real-time
PCR results was done on the data from H9 human embryonic stem cells to obtain the absolute expressions for
comparing more than two transcripts and comparative Ct method was done on the data from OVCAR8 ovarian
cancer cells and MCF7 breast cancer cells to obtain the ratio between a pair of transcripts.
H9 Stem cell line: Total RNA was extracted from human embryonic stem (ES) H9 cells by using TRIzol
(Invitrogen). To repeat the experiments of triplicate three times, 5µg RNA was used to synthesize complementary
DNA with ReverTra Ace (Toyobo) and oligo-dT (Takara) according to the manufacturer’s instructions. Transcript
levels of genes were determined by using Premix Ex Taq (Takara) and analysed with a CFX-96 Real Time system
(Bio-Rad). The templates for different transcripts were generated with PCR by using the template primers in
S1 Table in Supplementary. After isolation and purification, the templates were used to generate the standard
curves with qRT-PCR by using the qRT-PCR primers for different transcripts. The generated standard curves
have coefficient of determination (R2) over 0.999. The qRT-PCR primers were then applied to determine the
expression levels of different transcripts in H9 ES cells by calculating with the standard curves. The expressions
were carried out in three independent replications and the standard deviations were provided after the average.
Ovarian cancer cell line: 1µg of total RNAs were isolated from untreated OVCAR8 cells using Trizol (Invitrogen). RNA was reverse-transcribed using Superscript II reverse transcriptase (Invitrogen) according to
manufacture protocol. Real-time PCR was performed on CFX384 Real-time system (Bio-Rad) with FastStart
SYBR Green Master (Roche) with the primer sets in S2 Table in Supplementary. PCR conditions are 10 min
at 95◦ C and 40 cycles of 95◦ C for 45 sec and 60◦ C for 45 sec. Quantitation of the real-time PCR results was
done using comparative Ct method. Two replicates of qRT-PCR were performed using total RNAs isolated.
Breast cancer cell line: 0.5µg of total RNAs purified from MCF7 cells was used for oligo d(T)20 -primed
reverse transcription (Superscript III; Life Technologies). SYBR Green was used to detect and quantitate PCR
products in real-time reactions with the primer sets in S3 Table in Supplementary. PCR conditions for qRTPCR analysis are 2 min 94◦ C and 40 cycles of 94◦ C for 30 sec, 60◦ C for 20 sec and 72◦ C for 30 sec. Quantitation of
the real-time PCR results was done using comparative Ct method. GAPDH mRNA was used as a normalization
control for quantitation. Three replicates of qRT-PCR were performed using total RNAs isolated.

RNA-Seq data preparation
Three cell line RNA-Seq datasets were used for evaluating the accuracy of transcript quantification by comparison
with qRT-PCR results. The first dataset is the H9 embryonic stem cell line data from [28], downloaded from
SRA. The second dataset is an in-house dataset from the ovarian cancer cell line OVCAR8 prepared at University
of Kansas Medical Center. The third dataset is the MCF7 breast cancer cell line data from [29], downloaded
from SRA. There are 23,397,325 single-end 34bp reads in the stem cell line dataset, 19,892,473 paired-end 100bp

10
Cancer Type
Ovarian serous cystadenocarcinoma(OV)
Breast invasive carcinoma(BRCA)
Lung adenocarcinoma(LUAD)
Lung squamous cell carcinoma(LUSC)

Event
Survival
Relapse
Survival
Relapse
Survival
Survival

# of Patients by years
76(<3 ys) vs 62(>4 ys)
79(<1.5 ys) vs 68(>2 ys)
66(<5 ys) vs 57(>8 ys)
42(<5 ys) vs 38(>8 ys)
47(<2 ys) vs 56(>3 ys)
67(<2 ys)vs 77 (>3 ys)

Table 3. Summary of patient samples in TCGA datasets. The samples are classified by cutoffs on
survival and relapse time based on the available clinical information in each dataset.
reads in the OVCAR8, and 21,855,632 paired-end 76bp reads in the MCF7 mapped to the human hg19 reference
genome by TopHat2.0.9 [30] with up to 2 mismatches allowed. Exon coverages and read counts of exon-exon
junctions were generated by SAMtools [31] to be utilized with Net-RSTQ and base EM (equation (1)). Cufflinks
[32] directly infers transcript expressions based on the alignment by TopHat with the min isoform fraction set to
0 for better sensitivity.
TCGA RNA-Seq datasets of Ovarian serous cystadenocarcinoma (OV), Breast invasive carcinoma (BRCA),
Lung adenocarcinoma (LUAD) and Lung squamous cell carcinoma (LUSC) were analyzed for patient outcome
prediction with transcript expressions estimated by Net-RSTQ, base EM (equation (1)), RSEM [33] and Cufflinks
[32]. Both the gene expression and transcript expression data reported by RSEM [33] in TCGA (level 3 data)
were utilized as two baselines for cancer outcome prediction. The raw RNA-Seq fastq files (level 1 data) were
downloaded from Cancer Genomics Hub (CGHub) and processed by TopHat for use with Net-RSTQ, base EM
and Cufflinks. The patient samples in each dataset were classified into cases and controls based on the survival
and relapse outcomes as shown in Table 3. The command lines for preparing the data with RSEM and Cufflinks
are available in the S3 Text in Supplementary.

Results
There are six major results in this section, 1) isoform co-expression analysis on TCGA data to show the correlation
with protein domain-domain interactions; 2) overlapping the DDIs and KEGG pathways to understand the
transcript networks; 3) simulations for model validation and statistical analysis; 4) qRT-PCR experiments to
measure the performance of transcript quantification; 5) cancer outcome prediction on TCGA data to measure
the quality of transcript quantification as molecular markers; and 6) running time of Net-RSTQ.
Net-RSTQ was compared with base EM (the base model in equation (1)), Cufflinks [32] and RSEM (isoform
expression or gene expression) [33]. The accuracy of transcript quantification was directly measured on the
simulated data with ground-truth expressions and qRT-PCR data from the three cell lines. Cancer outcome
prediction on four TCGA cancer datasets evaluates the potential of using isoform expressions as predictive
biomarkers in clinical settings. Statistical assessment was also performed on randomized transcript networks to
evaluate the significance of the results.

Isoform co-expressions correlate with protein domain-domain interactions
To investigate the correlation between protein domain-domain interactions and isofrom transcript co-expressions,
we calculated the number of transcript pairs that are both nearby (being neighbors or having a distance up to
2) in the transcript network and highly co-expressed in the TCGA samples. The transcript co-expressions were
calculated by Pearson’s correlation coefficients of each pair of transcripts across all the samples in each dataset
with the isoform transcript quantification by Cufflinks. The transcript pairs were then sorted by the correlation
coefficients from the largest to the smallest and grouped into bins of size 1000. The number of transcript pairs

11

that are nearby in the transcript networks out of 1000 pairs are calculated within each bin and plotted in Fig 3(A)
and Fig 3(B) for the two cancer gene lists, respectively. In both Fig 3(A) and Fig 3(B), the left column shows
the plots of the number of pairs that are neighbors in the transcript network, and the right column shows the
plots of the number of transcript pairs with a distance up to 2 in the transcript network, among the 1000 pairs in
each bin. In all the plots, similar trends are observed in all the four cancer datasets: there are more interacting
isoform pairs in the bins with higher co-expressions. For example, among the 1000 transcript pairs with the
highest correlation coefficients, there are 73 interactions in the transcript network in OV dataset and thus, 73
interactions (y-axis) for bin index 1 (x-axis) is plotted in the left column of Fig 3(A). In all the plots, there is a
clear pattern that the numbers of matched nearby transcripts in the transcript network among the 1000 pairs in
the first few bins are higher than the expected average of 30 in the small network of density 3.02%, 114 in the
small network of density 11.41% (with distance up to 2), 45 in the larger network of density 4.54%, and 203 in the
larger network of density 20.33% (with distance up to 2). Moreover, the 2-step walk clearly promoted the number
of overlaps with the pairs of higher co-expressions in the small network. For example, the significant overlap is
extended from the first 25 bins to approximately the first 50 bins or more in the four datasets. The observation
suggests that higher co-expressions exist not only in the direct neighbors in the transcript network but also the
nearby nodes by a small distance. By exploring the network structure with prior information through neighbors
by many steps in iterations, Net-RSTQ model is expected to propagate the expression values from each transcript
to its nearby nodes in the network to capture the co-expressions. Note that considering the neighboring pairs
with distance up to 2 in the larger network will result in a graph of density 20.33%, which is likely to contain
too many false relations by the two-step walk. Thus, the plots of the larger network of distance-2 pairs are only
included for the completeness of the analysis.
The canonical 2x2 chi-square test was also applied to compare the number of the domain-domain interactions
in the first 10,000 transcript pairs (first 10 bins) with the number in the rest of the pairs. In all the four datasets in
both Fig 3(A) and Fig 3(B) with one exception in the LUSC dataset on the large network of distance-2 relation,
there is a significant difference that the highly co-expressed transcripts are more likely to interact with each
other in the transcript network, confirmed by the significant p-values. As explained previously, the exception
is likely due to the large number of false-positive pairs in the dense network. The observation further support
the hypothesis that protein domain-domain interactions correlate transcript co-expressions reported in previous
studies [12, 13].
To further understand the specificity of the domain-domain interactions in the highly co-expressed transcripts,
we calculated the number of domain-domain pairs that construct the DDIs in the top 10,000 co-expressed transcript pairs. The statistics suggest high diversity of the type of DDIs. For example, there are 547 interacting
transcript pairs among the 201 out of 898 transcripts in the top 10,000 co-expressed transcript pairs in OV dataset
for small network. The 547 interacting transcript pairs represent 770 different domain-domain interactions (There
might be more than one DDIs between a pair of transcripts). There are 739 interacting transcript pairs among
the 538 out of 5599 transcripts in the top 10,000 co-expressed transcript pairs in OV dataset for large network.
The 739 interacting transcript pairs represent 1277 different domain-domain interactions. The statistics suggest
that the correlation between protein domain-domain interactions and transcript co-expressions is not a bias due
to a few highly spurious DDIs. It is a general correlation in many different DDIs and co-expressed transcripts.
Very similar statistics were observed in all the datasets and both networks.
To further demonstrate the co-expression relations in the transcript network, two examples are shown in
S1 Fig. In S1(A) Fig, WHSC1L1 contains two isoforms connected with different interactions in the transcript
network. Isoform NM 017778 interacts with 12 transcripts with average correlation coefficients 0.22 and the
other isoform NM 023034 interacts with 13 more transcripts with average correlation coefficients 0.30 compared
with the average correlation coefficient 0.188 against the other unconnected isoforms across the samples in the
OV dataset. In S1(B) Fig, gene BRD4 contains two isoforms both of which are connected with the same 14
neighbors in the network. The average correlation coefficients between these two isoforms and the 14 neighboring
isoforms are both above 0.26 compared with the average correlation coefficient less than 0.15 against the other
unconnected isoforms across the samples on the BRCA dataset. In both examples, we observed high degree of

12

Figure 3. Correlation between transcript co-expression and protein domain-domain interaction
in TCGA datasets. The correlation coefficients between transcript expressions across all patient samples are
first calculated in each dataset for each pair of transcripts by Cufflinks. The correlation coefficients are then
sorted from largest to smallest and grouped into bins of size 1000 each. The x-axis is the index of the bins with
lower index indicating larger correlation coefficients. The y-axis is the number of the pairs among the 1000
pairs of transcripts in each bin coincide with protein domain-domain interaction between the transcript pair.
The red line is the smooth plot by fitting local linear regression method with weighted linear least squares
(LOWESS) to the curves. p-value is reported by chi-square test. (A) Co-expressions are calculated based on the
small gene list. (B) Co-expressions are calculated based on the large gene list. In both (A) and (B), the left
column shows the plots based on the connected transcript pairs in the transcript network and the right column
shows the plots based on the transcript pairs with distance up to 2 in the network.

13

agreement between co-expressions and DDIs.

Protein domain-domain interactions enrich KEGG pathways
To further understand the transcript networks, we overlapped the DDIs between genes in the two networks with
the 294 human KEGG pathways [34]. Among the 397 genes in the small network, 10.97%(17284) of the pairs are
co-members in at least one KEGG pathway. The 10.97% KEGG co-member pairs covers 42.70%(2122) of the
DDIs among the genes while the other 89.03%(140352) non-co-member pairs covers 57.30%(2748) of the DDIs.
By these numbers, there is about 6-fold enrichment of DDIs in the KEGG co-member genes in the small network.
Among the 2551 genes in the large network, the 5.15%(335372) KEGG co-member pairs covers 12.45%(40812)
of the DDIs among genes while the other 94.85%(6172229) non-co-member pairs covers 87.55%(287090) of the
DDIs. By these numbers, there is about 2.6-fold enrichment of DDIs in the KEGG co-member genes in the large
network. We also list the KEGG pathways that are highly enriched with DDIs in the large network in S4 Table.
Specifically, we consider the subnetwork of genes that are members of one KEGG pathway and calculated the
density of DDIs in the subnetwork to compare to the overall density of 5.04% in the whole network. Interestingly,
most of the enriched pathways are signaling pathways and disease pathways with very high DDI densities.

Net-RSTQ captures network prior in simulations
In the simulations, we applied flux-simulator [35] to generate paired-end short reads simulating real RNA-Seq
experiment in silico based on a ground truth transcript expression profile, using hg19 reference human genome
and RefSeq annotations downloaded from UCSC Genome Browser. To generate the ground-truth expression
profiles, the gene expressions were sampled from a poisson distribution and the proportions of the isoforms in
each gene were derived based on a neighbor average expression in the small transcript network and an initial mixed
power law expression profile with gaussian noise. A sequential updating was used to compute the proportion of
each isoform by adding the neighbors’ average expressions to the initial expression. The update procedure can be
found in the S2 Text in Supplementary. At last, flux-simulator was applied to simulate the short reads based
on the ground truth transcript expression file. 15 million 76-bp paired reads were generated by Flux Simulator
and mapped to the reference genome by TopHat [30] with up to two mismatches allowed. To account for the
large dynamic range of abundances, the expressions were normalized by log 2(expression+1).
The correlation coefficients between the transcript abundances estimated by Net-RSTQ under various λ, base
EM (equation (1)), Cufflinks and RSEM, and the ground truth transcript abundances are reported in Fig 4.
Furthermore, Net-RSTQ was also tested with 100 randomized networks with permuted indexes of transcripts
in the transcript network. To assess the impact of the network prior, two cases are shown. Fig 4(A) reports
the correlation between the transcripts in which isoforms coded by the same gene are connected with different
neighbors (109 out of 898 transcripts in 29 genes). Fig 4(B) reports the results from all the genes with more
than one isoform (712 out of 898 transcripts in 211 genes). In both comparisons, the transcript expressions
estimated by Net-RSTQ achieve higher correlation with the ground truth compared with base EM, Cufflinks and
RSEM. Slightly higher improvement was observed in the first case than in the second case since the network
prior plays more significant role in differentiating the isoform expressions by their different neighbors. When
randomized networks are used, Net-RSTQ leads to similar or worse results due to the wrong prior information.
Note that since the datasets were generated to partially conform to the network prior, the isoform expressions are
relatively “smooth” among the neighboring isoforms. Net-RSTQ tends to generate smoother expressions than
base EM, Cufflinks and RSEM. When applying Net-RSTQ with small λs and randomized network priors, slight
improvement was also observed due to the smoothness assumption on the data.

14

Figure 4. Correlation between estimated transcript expressions and ground truth in simulation.
In (A) and (B) x-axis are labeled by the compared methods and different λ parameters of Net-RSTQ. The bar
plots show the results of running Net-RSTQ with 100 randomized networks. In (C) and (D), x-axis are are the
percentage of edges that are removed from the networks. The plots show the results of running Net-RSTQ with
the incomplete networks. (A) and (C) report the results of 109 transcripts of the isoforms in the same gene
with different domain-domain interactions. (B) and (D) report the results of 712 isoforms in genes with
multiple isoforms.
To evaluate the effect of missing edges in the transcript network due to the undetected protein domain-domain
interactions, we randomly removed certain percentages of the edges in the transcript network and then run NetRSTQ with λ = 0.1 on the incomplete networks. The results are shown in Fig 4 (C) and (D) for the 109 transcripts
with different neighbors and the 712 transcripts in the gene with more than one transcript, respectively. It is
intriguing to observe that only when a large percentage of the edges are removed, the performance of NetRSTQ is affected. Intuitively, the observation can be explained by the fact that the Dirichlet prior parameter is
proportional to the average of the neighbors’ expressions. As long as some of the neighbors are still connected
to the target transcript in the network, the prior information is still useful. The result suggests that Net-RSTQ
is relatively robust to utilize transcript networks potentially constructed with a large percentage of undetected
protein domain-domain interactions.

Three qRT-PCR experiments confirmed overall improved transcript quantification
The isoform proportions estimated by Net-RSTQ, base EM, RSEM, and Cufflinks were compared to the qRTPCR results on the three cell lines. Parameter λ = 0.1 was fixed in all the Net-RSTQ experiments. Among the
genes that Net-RSTQ, base EM, RSEM, and Cufflinks report most different quantification results, qRT-PCR
experiments were performed to test the genes with relatively higher coverage of RNA-Seq data, coding two to
three isoforms, and the feasibility of designing isoform-specific primers in the qRT-PCR products (see S1, S2
and S3 Tables). Twenty-five genes in total were tested in the three cell lines: seven in H9 stem cell line, five in
OVCAR8 ovarian cancer cell line, and thirteen in MCF7 breast cancer cell line. The scatter plots of the relative
abundance of the first transcript in each gene estimated by Net-RSTQ, base EM, Cufflinks and RSEM were
compared to the qRT-PCR results in Fig 5(A) and (E). In the scatter plot, the estimated relative abundance

15

Figure 5. Validation by comparison with qRT-PCR results. (A) The scatter plots compare the
reported relative proportion of each pair of the isoforms of each gene between the computational methods
(Net-RSTQ, base EM, Cufflinks, and RSEM) and qRT-PCR experiments. The proportions of the two
compared isoforms in a pair are normalized to adding to 1. The x-axis and y-axis are the relative proportion of
one of the two isoform (the other is 1 minus the proportion) reported by qRT-PCR and the computational
methods, respectively. The scatter points aligning closer to the diagonal line indicate better estimations by a
computational method matching to the qRT-PCR results. The unshaded gradient around the diagonal line
shows the regions with scatter differences less than 0.1, 0.15, 0.2 and 0.25, within which the estimations are
more similar to the qRT-PCR results. (B)-(D) The scatter plots on each individual dataset. (E) The table
shows the percentage of predictions by each method within the unshaded regions and the overall Root Mean
Square Error of the predictions by each method compared to the qRT-PCR results.

16

by Net-RSTQ were closer to qRT-PCR results measured by the accuracy of various thresholds and Root Mean
Square Errors. Net-RSTQ achieved the lowest Root Mean Square Error of 0.291, which is more than 0.05 less
than 0.3435, the second best achieved by RSEM. In the 20% confidence region, Net-RSTQ puts 59.3% of the
pairs in the region compared with 37%, 29.6%, and 51.9% by base EM, Cufflink, and RSEM, respectively. RSEM
performed well by putting 37.0% of the pairs within 10% confidence regions but performed poorly in about half
of the pairs with more than 25% error.
The relative abundance of the seven genes in H9 stem cell line is shown in Fig 5(B), S2(A) Fig and S5 Table. In
all seven genes tested, the relative abundance estimated by Net-RSTQ is closer to the qRT-PCR results compare
to that by base EM and Cufflinks. RSEM performed similarly well on four genes and worse on the other three
genes, CBLC, TCF3 and NPM1. The same comparison on the five selected genes in OVCAR8 ovarian cancer cell
line is shown in Fig 5(C), S2(B) Fig and S6 Table. Cufflinks reports very low expressions in the first transcript
in four genes, three of which do not agree with the highly expressed transcript in the qRT-PCR results. While
base EM performed better for two genes (NSD1 and HNRNPA2B1), Net-RSTQ performed better on the other
three genes (HRAS, TSC2, and WHSC1L1). Net-RSTQ correctly predicted the overall enrichment of isoforms of
HNRNPA2B1 and NSD1 (NM 031243 > NM 002137 in HNRNPA2B1 and NM 022455 > NM 172349 in NSD1).
It is possible that the expressions of NM 002137 transcript in gene HNRNPA2B1 and NM 172349 in gene NSD1
were slightly over-smoothed by network information in Net-RSTQ with the fixed λ parameter. RSEM performed
slightly better on WHSC1L1 and NSD1 but much worse in the other three genes. The same comparison on the
thirteen genes in MCF7 breast cancer cell line is shown in Fig 5(D), S2(C) Fig and S7 Table. Cufflinks performed
poorly on 8 genes with more than 25% error while RSEM, base EM and Net-RSTQ performed poorly on 5, 4
and 3 genes, respectively. Overall, Net-RSTQ performed better than base EM and Cufflinks and slightly better
than RSEM. In summary, Net-RSTQ improved the overall isoform quantification significantly in the H9 stem cell
data and predicted more consistent cases in OVCAR8 and MCF7 cancer cell lines data. Note that there could be
more uncertainties in primer designs due to somatic DNA variations and cell differentiation and proliferation in
cancer cell lines, potentially a larger variation in the qRT-PCR experiments on the cancer cell lines is expected
than H9 stem cell line.

Net-RSTQ improved overall cancer outcome predictions
To provide an additional evaluation of the quality of transcript quantification, we designed six cancer outcome
prediction tasks by the assumption that better transcript quantification always leads to better isoform markers
for cancer outcome prediction. Net-RSTQ was compared with base EM, RSEM [33], and Cufflinks [32] by
classification with the quantification of isoform transcripts in two cancer gene lists (397 and 2551 genes) on
four cancer datasets. Each dataset is divided into four folds with two folds for training, one fold for validation
(parameter tuning), and one fold for test in a four-fold cross-validation. Support Vector Machine (SVM) with
RBF kernel [36] were chosen as the classifier. We repeated the four-fold cross-validation 100 times by each method
in each dataset.
The average area under the curve (AUC) of receiver operating characteristic of the 100 repeats are reported
in Table 4 when the small gene list was used and Table 5 when the large gene list was used. The transcript
expressions estimated by Net-RSTQ consistently achieved better average classification results than those by the
base EM. To evaluate the statistical significance of the differences between the AUCs generated by Net-RSTQ
and the base EM in the 100 repeats, we also report the p-values by a binomial test on the number of wins/loses
in all the experiments between Net-RSTQ and the base EM in Table 4 and Table 5. When the small gene list was
tested, three cases were significant with low p-values less than 0.001 and two cases were significant with p-values
just below 0.02 while in the BRCA (survival) data, the p-value is only moderately significant even though the
average by Net-RSTQ is higher. Overall, Net-RSTQ outperformed the base EM significantly. When the larger
gene list was tested, the improvements are not as significant. The improvement was only significant in one
dataset, BRCA (survival), and slightly significant in two datasets, OV (relapse) and LUSC (survival). In the
other three datasets, the improvements are not significant. Net-RSTQ also outperformed Cufflinks and RSEM
(transcript or gene) in five cases except the experiment on BRCA (relapse) dataset in Table 4. In Table 5, the

17

improvements are less obvious. Moreover, the isoform expression features are not more informative than gene
expression features. Overall, the classification performance with the small gene list in Table 4 is generally better
than or similar to the large gene list in Table 5 possibly suggesting less relevance to survival and relapse in the
large gene list.
Dataset
Net-RSTQ(Isoform)
base EM(Isoform)
RSEM(Isoform)
Cufflinks(Isoform)
RSEM(Gene)
p-value(Net-RSTQ vs base EM)

OV(Survival)
0.597
0.570
0.587
0.563
0.591
0.0011

OV(Relapse)
0.607
0.589
0.550
0.577
0.580
0.0198

BRCA(Survival)
0.683
0.673
0.651
0.676
0.651
0.1356

BRCA(Relapse)
0.590
0.542
0.616
0.593
0.558
2.248e-5

LUAD(Survival)
0.635
0.579
0.613
0.555
0.615
1.948e-8

LUSC(Survival)
0.567
0.550
0.536
0.556
0.559
0.0167

Table 4. Classification performance of estimated transcript expressions and gene expression on
the small cancer gene list. The mean AUC scores of classifying patients by estimated transcript (gene)
expression in four-fold cross-validation for each dataset are reported. The best AUCs across the five models
using isoforms as features are bold.

Dataset
Net-RSTQ(Isoform)
base EM(Isoform)
RSEM(Isoform)
Cufflinks(Isoform)
RSEM(Gene)
p-value(Net-RSTQ vs base EM)

OV(Survival)
0.599
0.590
0.584
0.562
0.604
0.3798

OV(Relapse)
0.585
0.572
0.569
0.582
0.577
0.0967

BRCA(Survival)
0.679
0.651
0.663
0.683
0.675
0.0018

BRCA(Relapse)
0.592
0.571
0.594
0.580
0.598
0.3822

LUAD(Survival)
0.604
0.597
0.587
0.583
0.627
0.6178

LUSC(Survival)
0.566
0.556
0.543
0.559
0.554
0.1356

Table 5. Classification performance of estimated transcript expressions and gene expression on
the large cancer gene list. The mean AUC scores of classifying patients by estimated transcript (gene)
expression in four-fold cross-validation for each dataset are reported. The best AUCs across the five models are
bold.
The parameter λ was tuned by the AUC on the validation set and the optimal λ was used to train the NetRSTQ model to be tested on the test set. The process is repeated for each fold in 100 repeats. To show the effect
of varying the λ on the classification performance in Net-RSTQ, we plotted the average AUC on the validation
set across the 100 repeats on the BRCA (survival) dataset with small gene list in S3(A) Fig. The optimal λ
was 0.1 in this experiment. The local gradient around the optimal λ suggesting that the transcript network is
playing an important role in inferring better transcript quantification from the RNA-Seq data. In S3(B) Fig, the
convergence of Net-RSTQ is also illustrated by each update through all the genes in each iteration. After less
than 10 overall iterations across 397 genes, Net-RSTQ converged well to a local optimum. Similar convergence
patterns were observed in all other TCGA samples.
To understand the role of the transcript network in the transcript expression estimation, we used 100 randomized networks to learn the transcript proportion in each experiment with λ fixed to be 0.1. In each randomization,
the edges were shuffled among all the transcripts in the small gene list. For transcript expressions learned by each
randomized network, we conducted the same four-fold cross validation to compute the average AUCs among 100
repeats. The boxplot of the AUCs learned with the 100 randomized networks is shown in Fig 6. Compared with
the classification results from the true transcript network, the result with randomized networks is always worse.
Another important observation is that, the median value of the AUCs across the 100 randomized networks is
lower or close to the result by the base EM, which suggests that the randomized networks play no role in improving classification and even lead to worse result. Overall, the results provide a clear evidence that the transcript
network is informative for the transcript expression estimation, and supplies more discriminative features for
cancer outcome prediction.

18

Figure 6. Statistical analysis with randomized networks. Comparison of the classification results by the
randomized networks and the true network. The λ parameter was fixed to be 0.1 in all the experiments. The
blue star and the red star represent the results with the real network and without network (base EM),
respectively. The boxplot shows the results with the randomized networks.

Running time
To measure the scalability of Net-RSTQ, we tested the Net-RSTQ algorithm on the data of the MCF7 breast
cancer cell line with three different networks, the small network (898 transcripts), the large network (5599
transcripts) and an artificial huge network (10000 transcripts). Fig 7 plots the CPU seconds of running NetRSTQ on the three networks under different λs. On the small network, the running time is at most about 100
seconds while on the large network and the huge network, the running time is in the scale of 1-e3 ∼1-e4 and
1-e5 ∼1-e6 , respectively. When λ = 0.1, the CPU time for the small network is 32.4 seconds; for the large network
is 2755 seconds; and for the artificial large network is 27806 seconds. The results suggest that Net-RSTQ might
scale up to about 10000 transcripts, and thus the performance is sufficient for studies focusing on any pathway
with up to several thousand genes in the pathway.

Figure 7. Running time. The plots show the CPU time (Intel Xeon E5-1620 with 3.70GHZ) for
running the Net-RSTQ algorithm one three networks, the small transcript network, the large
transcript network, and an artificial huge network of 10000 transcripts.

19

Discussion
In the paper, we explored the possibility of improving short-read alignment based transcript quantification with
relevant prior knowledge, protein domain-domain interactions. The observation of the correlation between isoform
co-expressions and protein domain-domain interactions suggests that the approach is a well-grounded exploration.
Different from previously methods [27], Net-RSTQ is a network-based approach that directly incorporates protein
domain-domain interaction information for transcript proportion estimation. The experiments suggested a great
potential of exploring protein domain-domain interactions to overcome the limitations of short-read alignments
and improve transcript quantification for better sample classification.
The Dirichlet prior from the neighboring isoforms play two different roles: differentiating isoform expressions
to reflect different functional roles or smoothing isoform expressions to reflect similar functional roles, depending
on whether the isoforms of a gene share the same or different interacting partners. This principle in modeling is
based on the hypothesis that isoforms playing different functional roles (e.g. containing different protein domains)
are more likely to behavior differently than isoforms with the same or similar functional roles (e.g. containing the
same protein domains). When the isoforms of a gene interact with different partners, their expressions correlates
with their partners’ expressions. And, when the isoforms of a gene interact with the same partners, there is no
benefit on differentiating their proportions to drive the functionality. A limitation is that when the functional
difference among the isoforms are not captured by domain content, the smoothing role might under-estimate
the difference in their proportions. Thus, our future goal is to bring in other type of functional information to
distinguish their functional roles in cancer such as preferential adoption of post-transcriptional regulations.
Currently, Net-RSTQ does not directly model multi-hits reads in multiple loci. In the TCGA experiments,
around 5-10% of the aligned reads in four datasets have multiple alignments reported by TopHat and only one of
the best alignments is considered. To check the effect of the multiple-alignment reads in transcript quantification,
we allow up to 20 best alignments by TopHat and normalized the read assignment qijk by the number of loci that
the reads aligned to. The correlation coefficients between the estimated gene expressions before and after the
normalization are above 0.98 in all the datasets. A potential rigorous solution is to add iteratively reassignment
of the reads to the potential origins based on updated abundance of the involved isoforms. The modification will
significantly decrease the computational efficiency and make it impractical on large RNA-Seq datasets.
There is also another alternative of integrating the network information directly as a regularization term on
the joint likelihood function of all the genes. We also explored this model in the S1 Text in Supplementary.
In the preliminary experiments, we observed very similar outputs between the alternative model and the NetRSTQ model shown in S8 Table. However, since the alternative model directly works with one large optimization
problem across all the genes, the convergence is much slower as shown in S4 Fig and the optimization package
used in the experiments ran into numerical issues. Thus, we believe the Net-RSTQ model is more scalable and
robust in comparison.
Currently, Net-RSTQ can scale on transcript network with up to around 5000 transcripts, which is sufficient
for more focused analysis of several thousand genes. The running time of Net-RSTQ on such large transcript
network is below 2 hours on each TCGA sample, compared with 5-8 hours needed for aligning the short reads.
To further scale up Net-RSTQ, we will investigate other faster strategies of utilizing short read information, such
as Sailfish [37] which directly estimates isoform expressions by counting k-mer occurrences in reads rather than
reads from the alignments. This will be our future direction.

Acknowledgments
The results are based upon data generated by The Cancer Genome Atlas established by the NCI and NHGRI.
Information about TCGA and the investigators and institutions who constitute the TCGA research network can
be found at http://cancergenome.nih.gov. The dbGaP accession number to the specific version of the TCGA
dataset is phs000178.v8.p7.

20

References
1. Mortazavi A, Williams BA, McCue K, Schaeffer L, Wold B (2008) Mapping and quantifying mammalian
transcriptomes by RNA-Seq. Nat Meth 5: 621-628.
2. Wang ET, Sandberg R, Luo S, Khrebtukova I, Zhang L, et al. (2008) Alternative isoform regulation in
human tissue transcriptomes. Nature 456: 470-476.
3. Li JJ, Jiang CR, Brown JB, Huang H, Bickel PJ (2011) Sparse linear modeling of next-generation mRNA
sequencing (RNA-Seq) data for isoform discovery and abundance estimation. Proceedings of the National
Academy of Sciences 108: 19867-19872.
4. Li W, Kang S, Liu CC, Zhang S, Shi Y, et al. (2014) High-resolution functional annotation of human
transcriptome: predicting isoform functions by a novel multiple instance-based label propagation method.
Nucleic Acids Research 42: e39.
5. Yang EW, Girke T, Jiang T (2013) Differential Gene Expression Analysis Using Coexpression and RNA-Seq
Data. Bioinformatics 29: 2153-2161.
6. Trapnell C, Hendrickson DG, Sauvageau M, Goff L, Rinn JL, et al. (2012) Differential analysis of gene
regulation at transcript resolution with RNA-Seq. Nature Biotechnology : 46–53.
7. Roberts A, Trapnell C, Donaghey J, Rinn J, Pachter L (2011) Improving RNA-Seq expression estimates
by correcting for fragment bias. Genome Biology 12: R22.
8. Li B, Ruotti V, Stewart RM, Thomson JA, Dewey CN (2010) RNA-Seq gene expression estimation with
read mapping uncertainty. Bioinformatics 26: 493-500.
9. Turro E, Su SY, Goncalves A, Coin L, Richardson S, et al. (2011) Haplotype and isoform specific expression
estimation using multi-mapping RNA-seq reads. Genome Biology 12: R13.
10. Huang Y, Hu Y, Jones CD, MacLeod JN, Chiang DY, et al. (2013) A Robust Method for Transcript
Quantification with RNA-Seq Data. Journal of Computational Biology 20: 167-187.
11. Jiang H, Wong WH (2009) Statistical inferences for isoform expression in RNA-Seq. Bioinformatics 25:
1026-1032.
12. Kriventseva EV, Koch I, Apweiler R, Vingron M, Bork P, et al. (2003) Increase of functional diversity by
alternative splicing. Trends in Genetics 19: 124 - 128.
13. Resch A, Xing Y, Modrek B, Gorlick M, Riley R, et al. (2004) Assessing the impact of alternative splicing
on domain interactions in the human proteome. Journal of Proteome Research 3: 76-83.
14. Tseng YT, Li W, Chen CH, Zhang S, Chen J, et al. (2015) IIIDB: a database for isoform-isoform interactions
and isoform network modules. BMC Genomics 16: S10.
15. Finn RD, Marshall M, Bateman A (2005) iPfam: visualization of protein–protein interactions in PDB at
domain and amino acid resolutions. Bioinformatics 21: 410-412.
16. Stein A, Ceol A, Aloy P (2011) 3did: identification and classification of domain-based interactions of known
three-dimensional structure. Nucleic Acids Research 39: 718-723.
17. Punta M, Coggill PC, Eberhardt RY, Mistry J, Tate JG, et al. (2012) The Pfam protein families database.
Nucleic Acids Research 40: 290-301.

21

18. Pruitt KD, Tatusova TA, Brown GR, Maglott DR (2012) NCBI Reference Sequences (RefSeq): current
status, new features and genome annotation policy. Nucleic Acids Research 40: 130-135.
19. Berman HM, Westbrook JD, Feng Z, Gilliland G, Bhat TN, et al. (2000) The Protein Data Bank. Nucleic
Acids Research 28: 235-242.
20. Yellaboina S, Tasneem A, Zaykin DV, Raghavachari B, Jothi R (2011) DOMINE: a comprehensive collection of known and predicted domain-domain interactions. Nucleic Acids Research 39: 730-735.
21. Prasad TSK, Goel R, Kandasamy K, Keerthikumar S, Kumar S, et al. (2009) Human Protein Reference
Database–2009 update. Nucleic Acids Research 37: 767-772.
22. Futreal PA, Coin L, Marshall M, Down T, Hubbard T, et al. (2004) A census of human cancer genes.
Nature Reviews Cancer 4: 177-183.
23. Higgins ME, Claremont M, Major JE, Sander C, Lash AE (2007) Cancergenes: a gene selection resource
for cancer genome projects. Nucleic Acids Research 35: D721-D726.
24. Xing Y, Yu T, Wu YN, Roy M, Kim J, et al. (2006) An expectation-maximization algorithm for probabilistic
reconstructions of full-length isoforms from splice graphs. Nucleic acids research 34: 3150-3160.
25. Li B, Ruotti V, Stewart RM, Thomson JA, Dewey CN (2010) RNA-Seq gene expression estimation with
read mapping uncertainty. Bioinformatics 26: 493-500.
26. Rozov R, Halperin E, Shamir R (2012) MGMR: leveraging RNA-Seq population data to optimize expression
estimation. BMC Bioinformatics 13: S2.
27. Pachter L (2011) Models for transcript quantification from RNA-Seq. ArXiv : 1104.3889.
28. Tanaka Y, Kim KY, Zhong M, Pan X, Weissman SM, et al. (2014) Transcriptional regulation in pluripotent
stem cells by methyl CpG-binding protein 2 (MeCP2). Human Molecular Genetics 23: 1045-1055.
29. Daemen A, Griffith O, Heiser L, Wang N, Enache O, et al. (2013) Modeling precision treatment of breast
cancer. Genome Biology 14: R110.
30. Kim D, Pertea G, Trapnell C, Pimentel H, Kelley R, et al. (2013) TopHat2: accurate alignment of transcriptomes in the presence of insertions, deletions and gene fusions. Genome Biology 14: R36.
31. Li H, Handsaker B, Wysoker A, Fennell T, Ruan J, et al. (2009) The Sequence Alignment/Map format
and SAMtools. Bioinformatics 25: 2078-2079.
32. Trapnell C, Williams BA, Pertea G, Mortazavi A, Kwan G, et al. (2010) Transcript assembly and quantification by RNA-Seq reveals unannotated transcripts and isoform switching during cell differentiation. Nat
Biotech 28: 511-515.
33. Li B, Dewey CN (2011) RSEM: accurate transcript quantification from RNA-Seq data with or without a
reference genome. BMC Bioinformatics 12: 323.
34. Kanehisa M, Goto S (2000) KEGG: Kyoto Encyclopedia of Genes and Genomes. Nucleic Acids Research
28: 27-30.
35. Griebel T, Zacher B, Ribeca P, Raineri E, Lacroix V, et al. (2012) Modelling and simulating generic
RNA-Seq experiments with the flux simulator. Nucleic Acids Research 40: 10073-10083.
36. Vapnik VN (1998) Statistical Learning Theory. Wiley-Interscience.

22

37. Patro R, Mount SM, Kingsford C (2014) Sailfish enables alignment-free isoform quantification from RNAseq reads using lightweight algorithms. Nat Biotech 32: 462–464.
38. Grant M, Boyd S (2014). CVX: Matlab software for disciplined convex programming, version 2.1. http:
//cvxr.com/cvx.
39. Grant M, Boyd S (2008) Graph implementations for nonsmooth convex programs. In: Blondel V, Boyd S,
Kimura H, editors, Recent Advances in Learning and Control, Springer-Verlag Limited, Lecture Notes in
Control and Information Sciences. pp. 95–110. http://stanford.edu/~boyd/graph_dcp.html.

23

Supporting Information
S1 Text. Alternative model by network-based regularization. We introduce a network-based regularization to the base model as an alternative model and evaluate the probabilities of a read being generated by the
transcripts in all the genes simultaneously as follows,
Lpen (P ; r) = log(L(P ; r)) − λkAP − W P k22 .

(16)

The term λkAP − W P k22 in equation (16) is a network constraint to encode prior knowledge from the transcript
network. Given a transcript interaction network, we assume that the connected transcripts are more likely to
co-express by introducing the following cost term over the expression π,

Ψ(P, π)

=

|T |
X

πi −

i=1

=

2


X
j∈nb(i)

πj 
|nb(i)|

2
X pj |rg(j) |
pi |rg(i) |


−
l
|nb(i)|l
i
j
i=1

|T |
X



|T |
X



j∈nb(i)

=

i=1

Aii pi −

2
X

Wij pj 

j∈nb(i)

= kAP − W P k22 ,

(17)

where nb(i) are the neighbors of transcript i. g(i) and g(j) are the genes containing transcripts i and j,
respectively. |rg(i) | denotes the number of reads aligned to gene g(i). li and lj are the length of transcripts
i and j. A is a diagonal matrix, where Aii = |rg(i) |/li . W contains the weights of transcript pairs in the
transcript network, where Wij = |rg(j) |/(|nb(i)|lj ). Minimizing Ψ(P, π) ensures that each transcript will receive
an expression close to the average expression of its neighbors in the transcript network. To solve equation (16)
we used CVX, a package for specifying and solving convex programs [38, 39]. The framework estimates the
expressions of transcripts in all the genes together in one optimization. We applied this framework to the small
network with 898 transcripts on MCF7 breast cancer cell line RNA-Seq data. Overall, the results between NetRSTQ and the alternative framework can be highly similar as shown in S8 Table when parameter λs are tuned.
However, the algorithm converges slowly (S4 Figure) compared to the convergence of Net-RSTQ in Figure 7 in
the main manuscript). It is clear that the alternative model does not scale to larger networks.
S2 Text. Steps of generating the simulation data. Below are the steps of generating the simulation data,
1. Generate 397 gene expression E by sampling a poisson distribution.
2. For the genes with multiple transcripts, generate isoform expressions from a mixed power law by the
flux-simulator to calculate the initial proportion pik (0) for each transcript Tik .
3. Let the initial isoform expression πik (0) = Ei ∗ pik (0) + gaussian noise
P
P
P
πg(j),j
4. Let α = 1; Repeat pik = (α ∗ j∈nb(Tik ) |nb(T
+ πik (0))/ q (α ∗ j∈nb(Tiq )
ik )|

πg(j),j
|nb(Tiq )|

5. After convergence, πik = Ei ∗ pik .
6. π is further normalized and used with flux-simulator as the ground truth expressions.

+ πiq (0)).

24

S3 Text. Cufflinks and RSEM command line. Cufflinks 2.1.1[32] was applied to generate isoform expression
as one of the baseline to compare with Net-RSTQ using the following command:
./cufflinks -p 4 -F 0 -G hg19RefSeq.gtf -o x.bam.
RSEM1.2.20[33] was applied to generate isoform expression as another baseline to compare with Net-RSTQ in
simulation and qRT-PCR studies using the following command:
./rsem-calculate-expression –paired-end –bowtie2 –bowtie2-path bowtie2/ -p 2 x1.fastq x2.fastq
hg19RefSeq x.

A

B

WHSC1
NSD1

NM_017778

SMARCA4

BRD4

WHSC1L1
0.22

MSH6

DNMT3A

MLL2

NM_014299

CREBBP

0.27

TRIM3
0.30

MLL3

PRDM1
EZH2

0.26

MLL
CBFA2T3

NM_023034

NM_058243

PBRM1 BRD3 EP300

SETD2

S1 Figure. Examples of transcript sub-networks with co-expression information. (A) Transcripts in
WHSC1L1 with correlation coefficients calculated on the OV dataset. (B) Transcripts in BRD4 with correlation
coefficients calculated on the BRCA dataset. Both examples are shown with the neighbors in the small transcript
network.

25

1

0.8
0.6
0.4
0.2
qRT−PCR Net−RSTQ base EM

RSEM

NM_004187

0.4

qRT−PCR Net−RSTQ base EM

NM_033360

0.6
0.4

qRT−PCR Net−RSTQ base EM

RSEM

0.4

TCF3

NM_001005862

1

RSEM

0.8
0.6
0.4

NM_023034

0.6
0.4

qRT−PCR Net−RSTQ base EM

RSEM

0

Cufflinks

qRT−PCR Net−RSTQ base EM

RSEM

Cufflinks

NM_080425

0.6
0.4
0.2

qRT−PCR Net−RSTQ base EM

RSEM

0

Cufflinks

NM_000267

qRT−PCR Net−RSTQ base EM

NM_022455

0.6
0.4

RSEM

Cufflinks

NSD1

1

NM_001128147

0.8

0

Cufflinks

0.8

NF1
1

RSEM

NM_016592

1

NM_172349

0.8
0.6
0.4
0.2

0.2

0

qRT−PCR Net−RSTQ base EM

GNAS
NM_004448

0.4

NM_017778

0.8

0.2

0.2

0

Cufflinks

0.6

0

Cufflinks

Relative abundance

Relative abundance

NM_001136139

RSEM

0.8

WHSC1L1

NM_003200

1

0.4
0.2

qRT−PCR Net−RSTQ base EM

0.2

qRT−PCR Net−RSTQ base EM

0.6

0.2

1

NM_175630

0.8

ERBB2

NM_004985

0.6

0

NM_153759

1

0.4

0.8

Cufflinks

NM_022845

0.6

0

Cufflinks

0.2

0.2

Relative abundance

RSEM

NM_001755

KRAS

1

NM_001146702
Relative abundance

Relative abundance

0.6

0

Cufflinks

0.8

0

0.8

DNMT3A

CBFB

0.8

0.2

KDM5C
1

NM_001130852

Relative abundance

0

NM_012116

Relative abundance

NM_005157

Relative abundance

NM_007313

Relative abundance

Relative abundance

1

C

CBLC

Relative abundance

ABL1

Relative abundance

A

qRT−PCR Net−RSTQ base EM

RSEM

0

Cufflinks

qRT−PCR Net−RSTQ base EM

RSEM

Cufflinks

NPM1
NM_001037738

PDGFB

NM_199185
1

0.8
0.6
0.4
0.2
0

RSEM

0.6
0.4

0.6
0.4
0.2

qRT−PCR Net−RSTQ base EM

RSEM

0

Cufflinks

qRT−PCR Net−RSTQ base EM

0.4
0.2

0.8
0.6
0.4

RSEM

0

Cufflinks

qRT−PCR Net−RSTQ base EM

RSEM

0.3
0.2

0.6
0.4

0.6
0.4
0.2

0

0

Cufflinks

Iso Group1

RSEM

RSEM

Cufflinks

Cufflinks

U2AF1
NM_001025203

0.4

NM_006758

0.8
0.6
0.4
0.2

qRT−PCR Net−RSTQ base EM

RSEM

0

Cufflinks

qRT−PCR Net−RSTQ base EM

RSEM

Cufflinks

WHSC1L1

NM_017778

0.8

0.6
0.4
0.2
qRT−PCR Net−RSTQ base EM

qRT−PCR Net−RSTQ base EM

1

0.6

0

0.8

0

0

Cufflinks

0.8

Relative abundance

Relative abundance

NM_023034

RSEM

Iso Group2

WHSC1L1
1

0.4
0.2

qRT−PCR Net−RSTQ base EM

0.2

qRT−PCR Net−RSTQ base EM

0.6

TP53
1

NM_001077183

0.8

0.2
RSEM

NM_000548

NM_003200

0.8

Cufflinks

Relative abundance

0.8

qRT−PCR Net−RSTQ base EM

0.4

TSC2
1

NM_172349
Relative abundance

1

0.5

0

NSD1
NM_022455

0.6

0.1

0.2
qRT−PCR Net−RSTQ base EM

NM_001136139

NM_020975
Relative abundance

0.6

NM_005343

NM_020630

Relative abundance

0.8

0

NM_176795

1

Cufflinks

1
0.7
Relative abundance

NM_002137
Relative abundance

NM_031243

RSEM

TCF3

RET

HRAS

HNRNPA2B1
1

NM_022768

0.8

Cufflinks

B
Relative abundance

NM_001201545

0.2

qRT−PCR Net−RSTQ base EM

RBM15

1
NM_033016

0.8

0

Relative abundance

NM_002608

Relative abundance

NM_002520

Relative abundance

Relative abundance

1

RSEM

Cufflinks

NM_017778

NM_023034

0.6
0.4
0.2
0

qRT−PCR Net−RSTQ base EM

RSEM

Cufflinks

S2 Figure. Evaluation by qRT-PCR experiments. The relative abundance of the transcripts in 7 tested
genes in H9 stem cell line (A), 5 tested genes in OVCAR8 ovarian cancer cell line (B), and 13 tested genes in
MCF7 breast cancer cell line (C) estimated by Net-RSTQ, base EM, Cufflinks and RSEM was compared with
the qRT-PCR experiments. The total abundance is normalized to 1 over the measured transcripts in each gene.

26
B

A

4

x 10
−5.14
−5.16

Log−Likelihood

AUC on validation Set

0.715

0.71

0.705

0.7

−5.18
−5.2

Log−Likelihood for each update
End of each iteration

−5.22
−5.24
−5.26
−5.28

10-4

10-3

10-2

0.1

0.2

0.3

0.4

0.5

−5.3
0

λ

500

1000

1500

2000

2500

Gene Updates

3000

S3 Figure. Model selection and convergence. The experiment done on BRCA (survival) dataset. (A)
Effect of varying λ on the classification performance. The plot shows the average AUC learned from the 100
repeats on validation set for different λs with the optimal λ in blue. (B) Convergence analysis by the total
log-likelihood. The plot shows the change of total log-likelihood in Net-RSTQ with each gene update. Each red
cross indicates the end of each round t in line 2 of Algorithm 1.

A

B

5

3.5

x 10

5

−3.5

CPU Seconds

Objective Function

x 10

−4

λ =1

−4.5
−5

2.5
2

−5.5
−6

3

0

500

1000

Iterations

1500

2000

1.5

0.1

1

10

λ

100

1000

S4 Figure. (A) Convergence and (B) Running time of the alternative regularized framework with
2000 iterations on MCF7 breast cancer cell line.
Gene/Transcript(isoform) Names
ABL1
NM 007313(iso2)
NM 005157(iso1)
CBCL
NM 012116(iso1)
NM 001130852(iso2)
KDM5C
NM 004187(iso1)
NM 001146702(iso2)
TCF3
NM 003200(iso1)
NM 001136139(iso2)
WHSC1L1
NM 023034(iso1)
NM 017778(iso2)
KRAS
NM 033360(iso2)
NM 004985(iso1)
NMP1
NM 002520(iso1)
NM 0010337738(iso3)
NM 199185(iso2)

Primer Names*
Template 1&2
qPCR 1&2
qPCR 2
Template 1&2
qPCR 1&2
qPCR 1
Template 1&2
qPCR 1&2
qPCR 1
Template 1&2
qPCR 1&2
qPCR 1
Template 1
Template 2
qPCR 1
qPCR 2
Template 1
Template 2
qPCR 1
qPCR 2
Template 1&2&3
qPCR 1&2&3
qPCR 1&3
qPCR 3

Forward
5-GGTTGGTGACTTCCACAGGAAA
5-TGAAAAGCTCCGGGTCTTAGG
5-TAGCCAAAGACCATCAGCGTT
5-ACCCTGTGGAACCAGGCTGC
5-ACCACCATTGACCTCACCTGC
5-CATCCTGCAGACCATCCCTG
5-GACCTGCTCGAGGTGACCCT
5-GCCTCTAACCAGCATTCCCA
5-AGAGGCTGAGGAGGTCCAGG
5-TGAATCCCAAAGCAGCCTG
5-TGAATCCCAAAGCAGCCTGT
5-GTATGCCTCCGTGGGACGA
5-CAGTTCCTCAGGCTACAGTGAAGA
5-CAGTTCCTCAGGCTACAGTGAAGA
5-GTCGGCGGCTTGATAAACAGT
5- CCCTTCAGCTACTGCAGATGC
5-TTCCTGCTCCATGCAGACTGT
5-TACATTGGTGAGGGAGATCCGA
5-TTCCTGCTCCATGCAGACTGT
5- TACATTGGTGAGGGAGATCCGA
5-TCCTTTCCCTGGTGTGATTCC
5-TCCTTTCCCTGGTGTGATTCC
5-AGCTGAAGAAAAAGCGCCAGT
5-AAGCCCAAAGATGGGGAGAA

Reverse
5-CACCGTCAGGCTGTATTTCTTCC
5-TTGACTGGCGTGATGTAGTTG
5-TTCGCGGTTATCAATTTCATGT
5-CACCTGCCCCAGCTCCAACT
5-ACTGCCAGGAGCTGCCAGTT
5-GGCCGAGCTCAGTCAGGTCT
5-AAGCTTTCTTCAGATCACAGGGAG
5-TCTCTGGAATGGTGATGGCC
5-CCAAGCCATTCTGGTTCTCC
5-TCTTGTAACTAATGTTTTTATTTTCCTTA
5-GGTTGTGGGCTTCGCTCAG
5-GGAGCTCCTGGACCCAGTGT
5-CATACAACAAACAGACATCTAGATCAAC
5-GTAATGTAGTTTCTTGCCAGCTTTACA
5-GTACCCATCCAGCTCAAACCG
5-CCAGGCACTCCAGGTGAAAGT
5-TAAGAAGTAATCAACTGCATGCACCA
5-TAAGAAGTAATCAACTGCATGCACCA
5-GCACCAAAAACCCCAAGACAG
5-TAGAAGGCATCATCAACACCCA
5-CATTGTCAGGTGAGGCAAATGC
5-TCGGCCTTTAGTTCACAACCG
5-CTTTTGTGCATTTTTGGCTGG
5-AAGGGCAAGGTTCACTGAATCA

S1 Table. Primer sets of the transcripts in seven genes of H9 stem cell line. * The numbers refer to
the isoforms in the first column.

27
Gene Name
HNRNPA2B1
HRAS∗
NSD-1
TSC2∗
WHSC1L1

Transcript Name
NM 031243
NM 002137
NM 176795
NM 005343
NM 022455
NM 172349
NM 000548
NM 001077183
NM 023034
NM 017778

Primer Sequence - Forward
TCC GCG ATG GAG GAA AAC TTT AG
AGC GGC AGT TCT CAC TAC AG
CCG CTC TGG CTC TAG CTC
AGG ATG CCT TCT ACA CGT TGG
TCG CCA TTC TTG CCA TTA GC
ATT GTC TGC TGC CCT TTT CC
CTC TCC ACC CGT GAA AGA ATT C
AAC GAG AGA CCC AAG AGG AT
ATG TAA AAC TGG GGC AGC AC
TTT CGG TTT GAG CTG GAT GG

Primer Sequence - Reverse
GCC ACC AAT AAA GAG CTT ACG G
TCC TTT TCT CTC CTC CAT CG
ACC AAC GTG TAG AAG GCA TCC
CAT GTC CTG AGC TTG TGC CT
TTT TCA TTG CTG CCG TCC AC
TGG AAT CTG GAT CAT CCC GA
GAC CAC ATG TTC AGA CAC ACT G
GA CGT ATC GAG CCA TCA TGT C
AAG CAC CAA CAG AAC AAC GC
TTT GGG CTG TTT GGC AAA CC

S2 Table. Primer sets of the transcripts in five genes of OVCAR8 cancer cell line. * Gene contains
more transcript(s) which can not be quantified by qRT-PCR.

Gene Name
ERBB2
NSD1
U2AF1∗
PDGFB
DNMT3A∗
GNAS∗
RBM15
RET
TCF3
WHSC1L1
CBFB
TP53
NF1∗

Transcript Name
NM 001005862
NM 004448
NM 022455
NM 172349
NM 001025203
NM 006758
NM 002608
NM 033016
NM 153759
NM 175630
NM 016592
NM 080425
NM 001201545
NM 022768
NM 020630
NM 020975
NM 001136139
NM 003200
NM 017778
NM 023034
NM 001755
NM 022845
Iso Group1
Iso Group2
NM 000267
NM 001128147

Primer Sequence - Forward
5-CACAGATAAAACGGGGGCAC
5-GAGGGCTGCTTGAGGAAGTAT
5-GACACGGTGCAGTCAAATCG
5-AGAAGAAATTGTCTGCTGCCC
5-TTGGAGCATGTCGTCATGGAG
5-TGCCCTCTTGAACATTTACCGT
5-CTCCGCGCTTTCCGATTTTG
5-GACTGAGCAGGAATGGTGAGAT
5-GCAGCTACTTCCAGAGCTTCA
5-GGGCAGCAGATACCCTGTTT
5-CGAGTCTTAGGCTGCGGAAT
5-CACTCCCGTCAACATGGACA
5-ATGCCTTCCCACCTTGTGAG
5-AACAAGAAGAGAGAAAACTTGGCG
5-TGCCCAGCAACTTAGGATGG
5-AATGGAAAGTCTACCGGCCC
5-GAGAAAGACCTGAGGGACCG
5-CAACTGCACCTCAACAGCGA
5-GCCTCTCAGTACAGCACTCC
5-AGAAAGGTGCCAGCGAGATT
5-GGATGCATTAGCACAACAGGC
5-CGGGAGGAAATGGAGGCAAG
5-GATGAAGCTCCCAGAATGCC
5-GAGGTGTAGACGCCAACTCT
5-TGAGGAAAACCAGCGGAACC
5-GTGGAATCCTGATGCTCCTGT

Primer Sequence - Reverse
5-CAGGGTCTGAGTCTCTGTGCT
5-TTTCTCCGGTCCCAATGGAG
5-GCTGCCGTCCACTTCATTTC
5-GGATCATCCGAAAGGGCTGT
5-CTGTGCACTGTTTTGGGGATT
5-CTGCATCTCCACATCGCTCA
5-AGAGGAAAAGGAACACGGCA
5-TCAAAGGAGCGGATCGAGTG
5-TTTCAGGCTACGATCCACGC
5-GGCTGGGCAGTACACAGAAT
5-GCACCTACCTTCCTGACCAC
5-GTACCCCGGAGAGGGTACTT
5-TCAACCAGTTTTGCACGGAC
5-TTTCCTCCCTTTAGGGACACC
5-TTGATTCCCACCCCAGAAGC
5-CAGAGCTCTTACCCGGTGTG
5-GGCCTCGTTAATATCCCGCA
5-CTCCAAGTTCAGGATGACCGA
5-GCCTGCCCATGTTAATGCTG
5-GCAGGTCACTCAGTCCTCTA
5-GCCAGCAGCTGTGAAACTCT
5-GTAAAGATGGGCAGCACACAT
5-GTAGCTGCCCTGGTAGGTTT
5-AAGTCAGGGCACAAGTGAACA
5-GCTGGCTAACCACCTGGTATAAA
5-AAAACCATAAAACCTTTGGAAGTGT

S3 Table. Primer sets of the transcripts in thirteen genes of MCF7 cancer cell line. * Gene contains
more transcript(s) which can not be quantified by qRT-PCR.

28
KEGG Pathway
Long-term potentiation
Fc epsilon RI signaling pathway
Fc gamma R-mediated phagocytosis
Axon guidance
VEGF signaling pathway
GnRH signaling pathway
Phosphatidylinositol signaling system
ErbB signaling pathway
Inositol phosphate metabolism
Neuroactive ligand-receptor interaction
Glycerophospholipid metabolism
Olfactory transduction
Neurotrophin signaling pathway
Central carbon metabolism in cancer
Natural killer cell mediated cytotoxicity
Glutamatergic synapse
Retrograde endocannabinoid signaling
RNA polymerase
Glioma
Ras signaling pathway
Inflammatory mediator regulation of TRP channels
Long-term depression
Gap junction
B cell receptor signaling pathway
Tight junction
Glycerolipid metabolism
mTOR signaling pathway
T cell receptor signaling pathway
Vascular smooth muscle contraction
Non-small cell lung cancer
Chemokine signaling pathway
Cell adhesion molecules (CAMs)
Choline metabolism in cancer
MAPK signaling pathway
Epithelial cell signaling in Helicobacter pylori infection
Oxytocin signaling pathway
Type II diabetes mellitus
Adherens junction
Lysine degradation
Circadian entrainment
Prolactin signaling pathway
Insulin signaling pathway
Cholinergic synapse
Morphine addiction
Bladder cancer
Platelet activation
Progesterone-mediated oocyte maturation
Toll-like receptor signaling pathway
Rap1 signaling pathway
Leukocyte transendothelial migration
Pancreatic cancer
Complement and coagulation cascades
Acute myeloid leukemia
Renal cell carcinoma
Endometrial cancer
Sphingolipid signaling pathway
Focal adhesion
Regulation of actin cytoskeleton
FoxO signaling pathway
Amyotrophic lateral sclerosis (ALS)
Osteoclast differentiation
Melanoma
Serotonergic synapse
Gastric acid secretion
Hepatitis C
Dopaminergic synapse
Calcium signaling pathway
Pathogenic Escherichia coli infection
Influenza A
Estrogen signaling pathway
Chronic myeloid leukemia
Aldosterone-regulated sodium reabsorption
Ovarian steroidogenesis
NOD-like receptor signaling pathway
cAMP signaling pathway
Colorectal cancer
Proteoglycans in cancer
Prostate cancer
RIG-I-like receptor signaling pathway
Signaling pathways regulating pluripotency of stem cells
cGMP-PKG signaling pathway
Thyroid hormone signaling pathway
TGF-beta signaling pathway
Endocytosis
ECM-receptor interaction
Adipocytokine signaling pathway
Hepatitis B
Ubiquitin mediated proteolysis
MicroRNAs in cancer
Pathways in cancer

# of gene in KEGG
67
70
91
127
61
92
82
87
63
275
95
407
120
67
134
116
103
32
65
228
99
60
89
72
138
59
60
104
121
56
189
142
101
257
68
159
48
73
52
97
72
140
113
93
38
131
88
106
211
118
66
69
57
66
52
120
207
215
134
51
131
71
114
75
133
131
180
55
175
100
73
39
51
57
200
62
204
89
70
142
167
119
80
213
87
70
146
137
297
398

# of overlapped gene
30
51
45
56
39
53
26
70
14
38
12
15
88
38
47
21
25
11
46
119
45
25
42
45
49
13
35
80
43
45
86
32
56
156
36
63
23
48
11
28
52
75
42
21
34
57
63
56
119
51
58
23
48
55
41
66
103
95
94
20
82
50
29
24
73
44
51
22
76
51
63
22
18
31
82
52
130
70
35
84
56
62
51
84
23
38
98
66
109
260

Interactions
552
1426
1002
1526
736
1356
318
2272
86
626
62
92
3128
580
848
168
238
46
800
5182
722
220
620
706
830
58
418
2178
622
662
2408
332
1012
7830
414
1264
166
722
36
232
792
1634
512
126
322
892
1074
842
3722
680
870
136
586
754
414
1068
2600
2210
2102
94
1544
566
190
130
1194
426
570
104
1172
524
796
96
64
186
1298
498
3082
884
220
1250
542
662
446
1204
90
242
1510
654
1484
5702

Density
61.33%
54.83%
49.48%
48.66%
48.39%
48.27%
47.04%
46.37%
43.88%
43.35%
43.06%
40.89%
40.39%
40.17%
38.39%
38.10%
38.08%
38.02%
37.81%
36.59%
35.65%
35.20%
35.15%
34.86%
34.57%
34.32%
34.12%
34.03%
33.64%
32.69%
32.56%
32.42%
32.27%
32.17%
31.94%
31.85%
31.38%
31.34%
29.75%
29.59%
29.29%
29.05%
29.02%
28.57%
27.85%
27.45%
27.06%
26.85%
26.28%
26.14%
25.86%
25.71%
25.43%
24.93%
24.63%
24.52%
24.51%
24.49%
23.79%
23.50%
22.96%
22.64%
22.59%
22.57%
22.41%
22.00%
21.91%
21.49%
20.29%
20.15%
20.06%
19.83%
19.75%
19.35%
19.30%
18.42%
18.24%
18.04%
17.96%
17.72%
17.28%
17.22%
17.15%
17.06%
17.01%
16.76%
15.72%
15.01%
12.49%
8.43%

S4 Table. Overlapped KEGG pathways with large transcript network. We consider the subnetwork of
genes that are members of one KEGG pathway and calculated the density of DDIs in the subnetwork.

29

Gene Name
ABL1
CBLC
KDM5C
KRAS
NPM1
TCF3
WHSC1L1

Transcript Name
NM 007313
NM 005157
NM 012116
NM 001130852
NM 004187
NM 001146702
NM 033360
NM 004985
NM 002520 (Iso1)
NM 199185 (Iso2)
NM 001037738 (Iso3)
NM 003200
NM 001136139
NM 023034
NM 017778

Net-RSTQ
64.46%
35.54%
73.12%
26.88%
80.52%
19.48%
51.36%
48.64%
34.92%
29.09%
35.99%
78.31%
21.69%
39.98%
60.02%

Estimated Proportion
base EM
Cufflinks
94.45%
16.48%
5.55%
83.52%
93.23%
87.59%
6.77%
12.41%
99.22%
99.95%
0.78%
0.05%
80.23%
58.82%
19.77%
41.18%
55.84%
0%
6.56%
53.97%
37.60%
46.03%
98.11%
96.42%
1.89%
3.58%
73.61%
4.50%
26.39%
95.50%

RSEM
53.14%
46.86%
87.54%
12.46%
91.52%
8.48%
19.67%
80.33%
84.62%
1.52%
13.86%
90.06%
9.94%
37.72%
62.28%

qRT-PCR
Results
56±4.4%
44±4.4%
51±9.8%
49±9.8%
86±4.5%
14±4.5%
36±4.2%
64±4.2%
52%∗
45%
3.2%
56±6.5%
44±6.5%
46±6.1%
54±6.1%

S5 Table. qRT-PCR results on H9 stem cell line. * Standard deviation of Iso1+Iso3 is 5.7% and Iso3 is
4.4%

Gene Name
HNRNPA2B1
HRAS∗
NSD1
TSC2∗
WHSC1L1

Transcript Name
NM 031243
NM 002137
NM 176795
NM 005343
NM 022455
NM 172349
NM 000548
NM 001077183
NM 023034
NM 017778

Net-RSTQ
60.07%
39.93%
54.92%
45.08%
55.40%
44.60%
52.39%
47.61%
48.84%
51.16%

Estimated Proportion
base EM
Cufflinks
81.42%
0%
18.58%
100%
100%
11.62%
0%
88.38%
96.39%
99.94%
3.61%
0.06%
92.36%
0.01%
7.64%
99.99%
63.85%
0.01%
36.15%
99.99%

RSEM
42.37%
57.63%
11.72%
88.28%
61.65%
38.35%
4.62%
95.38%
35.69%
64.31%

qRT-PCR
Results
82.83%
17.17%
50.16%
49.84%
99.98%
0.02%
45.36%
54.64%
18.81%
81.19%

S6 Table. qRT-PCR results on OVCAR8 cancer cell line. * Gene contains more transcript which can
not be quantified by qRT-PCR.

Gene Name
ERBB2
NSD1
U2AF1∗
PDGFB
DNMT3A∗
GNAS∗
RBM15
RET
TCF3
WHSC1L1
CBFB
TP53
NF1∗

Transcript Name
NM 001005862
NM 004448
NM 022455
NM 172349
NM 001025203
NM 006758
NM 002608
NM 033016
NM 153759
NM 175630
NM 016592
NM 080425
NM 001201545
NM 022768
NM 020630
NM 020975
NM 001136139
NM 003200
NM 017778
NM 023034
NM 001755
NM 022845
Iso Group1
Iso Group2
NM 000267
NM 001128147

Net-RSTQ
26.16%
73.84%
46.77%
53.23%
26.68%
73.32%
21.51%
78.49%
99.53%
0.47%
69.65%
30.35%
51.65%
48.35%
57.13%
42.87%
35.90%
64.10%
59.60%
40.40%
51.89%
48.11%
96.90%
3.10%
98.08%
1.92%

Estimated Proportion
base EM
Cufflinks
7.79%
3.35%
92.21%
96.65%
20.08%
21.94%
79.92%
78.06%
21.11%
73.15%
78.89%
26.85%
18.03%
43.93%
81.97%
56.07%
99.77%
99.45%
0.23%
0.55%
87.11%
93.62%
12.89%
6.38%
77.37%
63.21%
22.63%
36.79%
60.36%
70.78%
39.64%
29.22%
31.51%
6.70%
68.49%
93.30%
54.38%
77.61%
45.62%
22.39%
52.80%
13.84%
47.20%
86.16%
99.37%
99.10%
0.63%
0.90%
98.56%
15.14%
1.44%
84.86%

RSEM
6.34%
93.66%
17.14%
82.86%
25.33%
74.67%
18.12%
81.88%
99.46%
0.54%
89.66%
10.34%
73.85%
26.15%
67.76%
32.24%
36.31%
63.69%
62.43%
37.57%
57.67%
42.33%
97.41%
2.59%
94.55%
5.45%

qRT-PCR
Results
7.35±0.75%
92.65±6.0%
58.34±0.70%
41.66±0.75%
39.13±0.50%
60.87±2.5%
97.69±6.5%
2.31±1.4%
11.06±3.0%
88.94±4.0%
98.27±5.0%
1.73±0%
30.21±0.55%
69.79±1.8%
37.04±1.3%
62.96±2.3%
26.85±1.1%
73.15±4.0%
77.36±2.0%
22.64±0.18%
63.75±1.1%
36.25±1.2%
98.12±6.0%
1.88±0.30%
85.76±2.5%
14.24±2.0%

30

S7 Table. qRT-PCR results on MCF7 cancer cell line. * Gene contains more transcript(s) which can not
be quantified by qRT-PCR.

Net-RSTQ

λ
1e-04
1e-03
5e-03
0.01
0.05
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
0.99
0.999

Regularized framework
0.1
1
10
0.9689
0.9536
0.9442
0.9696
0.9549
0.9453
0.9722
0.9583
0.9488
0.9751
0.9620
0.9526
0.9887
0.9804
0.9720
0.9943
0.9895
0.9819
0.9959 0.9947
0.9882
0.9950 0.9954 0.9899
0.9938
0.9951
0.9905
0.9925
0.9945
0.9909
0.9910
0.9935
0.9912
0.9889
0.9920 0.9914
0.9854
0.9891
0.9911
0.9770
0.9815
0.9883
0.9328
0.9387
0.9543
0.9121
0.9181
0.9349

100
0.8668
0.8680
0.8715
0.8751
0.8939
0.9038
0.9111
0.9141
0.9164
0.9188
0.9218
0.9258
0.9319
0.9424
0.9601
0.9578

S8 Table. Correlation Coefficients between the results of Net-RSTQ and the alternative regularized
framework with different λs. The highest correlation coefficients for each λ in the alternative regularized
framework is bold.

