The data of CNA, RNA-seq data, and clinical data of CRC were downloaded from the TCGA database. By mapping the copy number probe across the reference genome of hg38, the SCNA at gene level was calculated using Gistic2 software (Mermel et al., 2011). The value of SCNA represents the portability of copy number alteration and the q-value for the genes in aberrant regions. The q-value > 0.1 and q-value < −0.1 were considered as copy number amplified and deleted, respectively. For each gene, the samples with SCNA value > = x (x represents the threshold of SCNA with a value >0) were identified as copy number amplification samples (CNAS), the samples with SCNA < = −x were identified as copy number deleted samples (CNDS), and the samples with | SCNA | < x were identified as copy number non-altered samples (CNNS). The location information of chromosomes was obtained from the HGNC database (Braschi et al., 2019). RNAseq FPKM data was downloaded from University of California Santa Cruz (UCSC, http://genome.ucsc.edu/), and more than 80% of genes with 0 value were filtered out. The test data-set was collected from the Cancer Cell Line Encyclopedia (CCLE; http://www.broadinstitute.org/ccle/home).

PSGs of SCNA were screened in five steps as described below:

In order to verify the stability of the dosage-sensitivity of PDSGs, the correlation coefficients between gene expression and copy number alteration were calculated with the RNA-seq of CRC and CNA at gene level downloaded from the CCLE database. These values were compared with the findings obtained from TCGA.

In order to further identify the genes affected by PDSGs, Pearson correlation coefficients of these six PDSGs and other genes was calculated as co-expression values in CNAS or CNDS, CNNS. Gene pairs with correlation coefficients higher than 0.5 in one group and less than 0.1 in another group were screened as differentially co-expressing gene pairs. Network visualization tools were executed using Cytoscape (Shannon et al., 2003).

All the analysis was performed in the R computing environment. Survival curves were estimated using the Kaplan-Meier method. Gene function enrichment was performed using the Cluster Profiler package (Yu et al., 2012).

A total of 448 CRC samples with SCNA and RNA-seq data were downloaded from The Cancer Genome Atlas (TCGA). The samples were screened for survival information. There were 22,752 genes, of these 17,442 were protein-coding and 14,688 were differentially expressed.

After applying FDR < 0.1 and FC > 1.2, 6,814 genes had up-regulated expression in CNAS. Twenty-five genes had a down-regulated expression in CNDS. Cox regression analysis was applied to calculate the correlation between SCNA and survival time. A total of 215 prognosis-sensitive genes (PSGs) significantly related to SCNA were obtained, of these 214 were amplified and one was deleted. Next, the 21 SCNA threshold value was raised from 0.1 to 0.5 at a step of 0.02. For each threshold, the samples were classified into CNNS, CNAS, CNDS group and logRank test between CNNS and CNAS, CNDS and CNNS was performed. As shown in Figure 1 , 73.02% of genes didn’t show any significant classification with any threshold. A total of 15 genes showed stable prognosis classification of patients in more than 10 threshold values, suggesting these 15 genes can be considered as stable markers for prognosis classification in CRC.

After further screening stable PSGs which are highly affected by copy number dosage effect, the Pearson correlation coefficient between copy number and corresponding expression value (FPKM) of these 15 genes was calculated. Finally, six genes (NDUFB4, WDR5B, IQCB1, KPNA1, and SEC22A) which are stable PSGs ( Figure 2 ) were identified. The average dosage effect score was 0.5918 and the variance was 0.066.

Kaplan-Meier survival curve analysis revealed six (6) PDSGs with similar results in a different threshold of SCNA. In the 0.1 SCNA threshold value, genes GTF2E1, NDUFB4, IQCB1, KPNA 1 and WDR5B had a significant classification effect ( Figures 3A–C ). At the 0.3 threshold value of SCNA, all six genes had a similar and significant classification effect ( Figure 3D ). At the 0.5 threshold value, five genes (GTF2E1, NDUFB4, IQCB1, KPNA1, WDR5B) had similar classification effect ( Figures 3E, F ). Although the statistical significance of the two classifications (p-value = 0.087199 and p-value = 0.12643) in 0.5 SCNA threshold was not significant, their classification curves were distinctly separated. The non-significance can be primarily attributed to the very small number of samples with SCNA threshold >0.5.

In order to verify if the copy number of six PDSGs is dosage-sensitive in the data from cell lines with 53 cell line samples, the dosage effect score of these six PDSGs in CRC from CCLE was calculated. An average score of 0.5978 and variance was 0.082 consistent with the result from TCGA was obtained ( Figure 4A ). The Pearson correlation coefficient was 1, suggesting that the gene dosage effect is stable in CRC different data.

Further to test similarity between survival curves of these six PDSGs, we mapped them to chromosomes and found that they all are located on 3q13.33–3q21.1. By computing the correlation coefficients between the copy number of two pairs of genes an average value of 0.9967 ( Figure 4B ) was observed. This indicates that these six PDSGs are highly consistent with each other during alteration.

Research have shown that heterogeneity of copy number alterations exists in ongoing unstable chromosome in COAD (Bolhaqueiro et al., 2019). There are some chromosomes fragile sites in genome, the genes in fragile sites may break when they fell external pressure. In order to determine the presence of breakpoints in the region near to 6PDSGs, they were mapped on the database of human chromosomes fragile sites (HumCFS, http://webs.iiitd.edu.in/raghava/humcfs/). As a result, FRA3D (3q25.32) and FRA3C were found to be near to six PDSGs. Correlation analysis of SCNA in six PDSGs and the genes in FRA3D and FRA3C was performed. Gene RSRC1 (R = 0.82), MLF1(R = 0.82) in FRA4D, and LPP (R = 0.80) in FRA3C had lowest relationship with PDSGs. Thus we infer that the breakpoints in fragile site may explain the reason for the nearby region and a similar SCNA value.

In order to further explore if these six PDSGs can also affect the expression of other genes in CRC, we screened genes with (R) > 0.5 and (R) < 0.1 in a different class of samples by calculating the differences of gene co-expression between CNAS and CNNS. A total of 234 co-expressed gene pairs were observed and 215 genes ( Figure 5A ) involved in differential co-expression networks were identified. The whole network constitutes a component suggesting that CRC is a disease involving multiple genes. Among these 194 gene pairs were co-expressed in alteration samples (R > 0.5), but not co-expressed in non-alteration samples (R < 0.1), while the other 40 pairs behaved in a reverse manner. In the network, gene NDUFB4, SEC22A had the highest degree (109 and 45 respectively) consisting of 15 co-linked genes. The genes CAPN14 and CMPK2 were affected by three PDSGs (NDUFB4, SEC22A, and IQCB1). This suggests that PDSGs are closely linked and interact with each other.

Each PDSG in the network was related to at least 13 genes and 22 genes were associated with more than one PDSG. We also found that several PDSGs-associated genes were also COAD-related. The co-expression of GTF2E1-WNT8B was activated in CNAS(R = 0.59). WNT8B one member of the WNT signal was differentially expressed in COAD (Neumann et al., 2014). In addition to this, after mapping the PDSG-related genes to the driver gene list from DriverDB (Liu et al., 2019), three genes (C8orf33, LAPTM4B, PTP4A3) were found ( Figure 5B ), and they all were co-expressed with gene NDUFB4 in CNAS but not in CNNS. Mapping of PDSG-related genes on the tumor suppressor database (TSGene, http://bioinfo.mc.vanderbilt.edu/TSGene/) revealed 16 TSGs ( Figure 5A , Triangle). Among these, gene DCDC2, ISG15, RARRES3 can affect more than one PDSG. Gene RARRES3 has been shown to be mutated, differentially expressed and also inhibits metastasis in COAD (Lee et al., 2018a). ISG15 is shown to have significant differential expression in COAD (Yu et al., 2019; Zamanian-Azodi and Rezaei-Tavirani, 2019).

Further to explore the possible functions of these six PDSGs, linked genes were extracted and gene ontology function enrichment analysis was performed. Genes linked to gene NDUFB4 ( Figure 5C ) were mainly enriched in functions such as “transmembrane receptor,” “transmembrane transport,” “peptide receptor,” “G protein-coupled receptor,” “transforming growth factor.” Genes linked to gene GTF2E1 were enriched ( Figure 5D ) in functions such as “cyclin-dependent protease,” “ATP synthase transport proton-related functions.” Previous studies have shown that transforming growth factor can also promote tumorigenesis (De Miranda et al., 2015; Yu et al., 2018; Kim et al., 2019). G-protein-coupled receptors (GPCRs) are a member of the largest cell surface molecule family involved in signal transduction and are considered as the key molecule in the growth and metastasis of tumors (Wielenga et al., 2015; Insel et al., 2018). Malignant cells often hijack the normal physiological functions of GPCRs to survive, proliferate independently, escape the epidemic system, increase blood supply, invading the surrounding tissues and spread to other organs.