使用说明：

Rscript $scriptdir/lasso_cox.r -h 
usage: /work/my_stad_immu/scripts/lasso_cox.r [-h] -i data -t time -e event -v
                                              variate [variate ...] [-s seed]
                                              [-l lambda] [-P predict.time]
                                              [-o outdir] [-p prefix]

lasso cox regression analysis

optional arguments:
  -h, --help            show this help message and exit
  -i data, --data data  input data file path[required]
  -t time, --time time  set suvival time column name [required]
  -e event, --event event
                        set event column name [required]
  -v variate [variate ...], --variate variate [variate ...]
                        variate for cox analysis [required]
  -s seed, --seed seed  set random seed [default 2021]
  -l lambda, --lambda lambda
                        set lambda cut off [default NULL]
  -P predict.time, --predict.time predict.time
                        Time point of the ROC curve to select cutoff [default
                        365 1095 1825]
  -o outdir, --outdir outdir
                        output file directory [default cwd]
  -p prefix, --prefix prefix
                        out file name prefix [default lasso_cox]

使用举例：

Rscript $scriptdir/lasso_cox.r -i imm.unicox.metadata-exp.tsv -e EVENT -t TIME \
  -v SYT12 CDH2 GPNMB TMIGD3 LINC01094 SLC22A20P IGHV4-61 IGHV2-5 SERPINA5 MS4A4A FAM83A IGLV3-9 STARD3    -o lasso

参数说明：

-i 输入生存数据与基因表达文件 

结果展示：

lasso变量筛选与模型构建：

 Lasso(Least absolute shrinkage and selection operator)方法进一步筛选变量优化模型。  Lasso主要用来进行变量筛选，特别是当自变量之间有共线性时很有用， 我们用10fold 交叉验证的方法排除共线性严重的基因优化和简化模型。该分析可以使用R包glmnet(Friedman et al. 2010)完成。

注：LASSO为寻找最佳的模型，引入变量λ (lambda 又叫收缩算子、模型系数比、调优系数或惩罚值)  随着λ增加，各变量的回归系数β在减小，有些会变为0，说明该变量在此时对模型贡献微乎其微，可以剔除。图中，一条彩线代表一个变量的回归系数β值的变化，x轴下方的数字为惩罚值（调优系数），x轴上方为在该值下的剩余的变量个数。LASSO 回归就是通过生成一个惩罚函数对回归模型中的变量回归系数进行压缩，达到防止过度拟合，解决严重共线性的问题。

λ值决定了哪些变量可以使模型最优，使用交叉验证可寻找最佳λ值。Partial-likelihood deviance (偏似然偏差) 随Log(λ)变化曲线，此值越小说明模型拟合越好。图中给出了两个惩罚值（调优系数）λ：  一个是当偏似然偏差最小时的λ值，即lambda.min；

We applied the Cox regression model with LASSO based on the R package “glmnet” to construct an optimal  gene‐associated prognostic model。

The Risk score was calculated with the following formula: The  risk score=

, where Expri represents the expression level of gene i and coefi represents the regression coefficient of gene i in the signature.We grouped all patients into low- or high-risk groups according to the median value of XXXX‐based risk signature and performed survival analysis with Kaplan-Meier method. The logrank test was used to compare the difference in the survival status between the high‐ and low‐risk groups.

风险评分

根据模型计算各样本分风险值，按照风险值的中位数将样本划分为高低风险组，分别绘制风险值分布散点图，生存时间散点图，signature基因表达热图。

高低风险组预后差异分析：绘制Kaplan-Meier生存曲线，并用Log Rank法检验两组的生存率是否有差异。

模型预测性能评估

模型的好坏可以从区分度（Discrimination）和一致性（Calibration）两方面考虑。区分度主要用于反映预测模型的区分能力，是评估模型有多大把握确定它所预测的患者发生该事件的能力。一致性指结局实际发生的概率和预测的概率的一致性或者接近程度。前者可通过ROC曲线下面积（AUC）或C统计量来评价，后者可通过校准图来评价。以下为模型ROC曲线：

To reflect the prediction ability of the XXXX‐based risk signature, we generated the time-dependent receiver operating characteristic curve (ROC) and calculated the area under the curve (AUC)  (R package “survivalROC” ) for 1-year, 3-year, and 5-year overall survival (OS). The Kaplan-Meier, log‐rank, ROC curve, and calibration analyses were all performed and visualized by the “survivalROC”, “rms”, “survival”, and “survminer” packages.

脚本获取与使用课程：https://study.163.com/course/introduction/1211864801.htm?share=1&shareId=1030291076

References