title: "Analysis of Liver Hepatocellular Carcinoma Data in R"

author: 'Bo Yang, Rui Nie, Yueying Hu'

output:

pdf_document:

fig_width: 3

fig_height: 2

latex_engine: xelatex

header-includes:

- \usepackage{geometry}

- \geometry{left=0.6in, right=0.6in, top=0.6in, bottom=0.6in}

- \usepackage{setspace}

- \singlespacing

- \usepackage{enumitem}

- \setlist[itemize]{noitemsep, topsep=0pt}

- \setlist[enumerate]{noitemsep, topsep=0pt}

fontsize: 8pt

bibliography: references.bib

```{r, setup, message=FALSE, include=FALSE}

```{r, message=FALSE, include=FALSE}

Hepatocellular carcinoma (LIHC/HCC) is the one of the leading cause of cancer death, which accounts for 8.2% of global cancer incidence and 4.7% of cancer-related mortality [@fitzmaurice2017global] Previous study suggested that lifestyles and personal habits, including drinking alcohol, smoking, are positively associated with the onset LIHC [@kuper2000tobacco], which has greatly enhanced the efficiency of screening for the disease by health providers. Nevertheless, the difficulty for disease detection remains alarmingly high. Therefore, an indicator with higher sensitivity would advance the diagnosis stage and prolong patients' survival.

RNA Sequencing (RNA-seq) using the capabilities of high-throughput sequencing methods provides insight into the transcriptome of a cell. Through understanding the gene expression patterns reflected by RNA sequencing results, we aim to locate indicator genes for LIHC with the help of statistical learning methods Furthermore, given the massive amount of existing genes in human genome, our problem can be formulated as a high-dimensional problem. One of our side goal is to perform variable selection to reduce the risk of overfitting.

We're primarily interested in analyzing RNA-seq data provided through The Cancer Genome Atlas (TCGA) for liver cancer patients. The dataset includes data from 424 LIHC participants and over 20,000 genes for expression level analysis. The inclusion of both normal and tumor samples allows us for comparative analysis between cancerous and non-cancerous liver tissues. We used `GDCquery`, `GDCdownload`, and `GDCprepare` functions for data accessing to download and prepare the LIHC data from TCGAbiolinks.

```{r, message=FALSE, results='hide'}

Normalization for reads count help account for gene length and sequencing depth. Fragments Per Kilobase of transcript per Million mapped reads (FPKM) provides a comparable measure to interpret and utilize RNA-seq results. We standardized our data into FPKM format for fair representation and comparison gene expressions.

We collected the sample IDs from normal samples and from tumor samples separately and retrieved their index for future train-test dataset split. In summary, out of 424 samples, there are 50 normal samples and 374 tumor samples. This suggests we consider additional metrics such as precision, F1, or AUC-ROC scores to account for potential bias introduced by the imbalanced labels aside from crude accuracy for classification tasks.

Of all the variables (genes) selected for our preprocessed data, we'd like to filter out the ones with constant values across samples, as they provide no information in differentiating tumor samples from normal samples. This leads to the removal of 381 genes from the entire sample.

To select variables while maintaining interpretation along the process, we will use the variance filtering method to select genes with a more variable expression activity from the samples we included. We kept the variance cutoff to be as high as 0.95, leaving only 5% (985) of the original genes for future statistical analysis.

For building predictive statistical learning methods in LIHC classification, we continue with 985 selected genes using Elastic Net, decision tree, and random forest and assess their performances under the idea of train-test split validation. We will also use cross-validation to tune the hyperparameters for each model and compare their performances in terms of specificity, sensitivity, precision, F1, and AUC-ROC scores. \## Train-Test Split Considering the imbalanced labels of two classes, it is important to preserve the representation of both classes in the training and testing samples to ensure the model is not overly biasing towards one class over the other. As a consequence, we use `createDataPartition` to create a stratified split for training and testing purposes with a 3:1 ratio.

Elastic Net is a penalized regression method that combines the L1 and L2 penalties of Lasso and Ridge regression. The L1 penalty suppress the coefficients of some variables to zero, which is equivalent to variable selection. The L2 penalty shrinks the coefficients of some variables towards a smaller scale. We take the idea of logistic regression and use below formula as the loss function for Elastic Net in classification problems:

We determine the alpha and lambda coefficient through grid search, where alpha determines the balance between L1 and L2 penalty and lambda determines the strength of the penalty.

```{r, warning=FALSE, message=FALSE}

Decision Tree is a non-parametric supervised learning method that can be used for both classification and regression problems. It works by recursively partitioning the data into smaller subsets based on the most significant variable at each step. The partitioning process is repeated until the data is completely separated or the stopping criteria is met. We displayed the gene names of the most important variables selected by Decision Tree.

Random Forest is an ensemble learning method that combines multiple decision trees to improve the performance of the model. It works by building multiple decision trees on different subsets of the data and then averaging the predictions of all the trees. The random forest model is more robust to overfitting and less sensitive to outliers than a single decision tree. Here, we use the `randomForest` function to build the random forest model and tune the hyperparameters by cross validation.

Further, we explore the the variable importance, we consider two types of measure of variable importance: Mean Decrease Accuracy and Mean Decrease Gini. We see the gene TSPAN6 is the most important predictor in Mean Decrease Accuracy and in Mean Decrease Gini, which is consistent with the results of the decision tree and elastic net regression.

Overall, the three statistical learning methods have similar performance in terms of specificity and raw accuracy, and they all successfully identify similar genes as the most important predictors ("TSPAN6", "TNMD", "DPM1", "SCYL3", "C1orf112", "FGR"). However, the single decision tree model has the lowest sensitivity, F1 score, and AUC-ROC score, which indicates that the model is overfitting. Comparably, the elastic net regression model and the random forest model are more robust to such unbalanced dataset. All code is avaliable in github: <https://github.com/rrrrn/LIHC625/tree/main>

```{r, warning=FALSE, message=FALSE, echo=FALSE, results='hide', out.width="50%"}

| Member | Contribution |

|----------|-------------------------------------------------------------|

| Bo Yang | Data Preprocessing, Data analysis, write LIHC data from TCGA and Preprocessing part in the report |

| Rui Nie | Variable Selection ,Data analysis, write Introduction, Variable Selection and Elastic Net part in report |

| Yueying Hu | Data collection, Data analysis, write Discussion, Decision Tree, Random Forest part in report |