Hello Biostars community,

I'm trying to analyze miRNASeq data of TCGA melanoma (SKCM) samples. This is a prototype analysis in which I picked `hsa-mir-155` and looked at how its expression is correlated with survival. My end goal is expanding this analysis to other cancer types also focusing on other miRNAs.

I performed my analyses with two of the popular packages in R (RTCGA and TCGAbiolinks) and obtained quite different results. I'm trying to make sense of what might be causing the different results. Any help is appreciated. For comparison purposes, I'm attaching some figures and my codes as well.

![Kaplan-Meier curves][1] Overall Kaplan-Meier curves look quite different. These plots were generated on the whole cohort (unsegregated for any parameter). I'm including the plot from `oncoLNC.org` database as well for comparison purposes. RTCGA looks more like the `oncoLNC` curve.


You can see similar differences when the data is segregated based on `patient.gender`:
![KM_curves_per_gender][2]


When I segregate patients based on the expression of `hsa-mir-155` (top and bottom thirds), the differences become more obvious:
![mir155_top_bottom_thirds][3]


To understand what might be different in the datasets I exported the data from both packages and performed a comparison. [Linked][4] excel file shows the comparisons including clinical details (`days_to_death` and `days_to_last_follow_up`) and gene expression values of `hsa-mir-155` (both `read_count` and `reads_per_million_miRNA`.  I noticed that there are considerable differences between two datasets. I'm pretty sure, the way I organized the data is ok and I don't think the differences are due a mistake in data manipulation in R.

The code I used for analyses can be found at
[RTCGA_v1][5] and 
[TCGABiolinks_v1][6]


Please let me know why you think there is a discrepancy here. I'm pretty new to this type of analyses and hopefully, I didn't miss something silly.

Thank you very much in advance for your insights,

Atakan


  [1]: https://s22.postimg.cc/f3513t4wx/Untitled.jpg
  [2]: https://s22.postimg.cc/h8zbxuq8h/untitled2.jpg
  [3]: https://s22.postimg.cc/ra4678vmp/Untitled3.jpg
  [4]: https://www.dropbox.com/s/acy1ubbk1fkdhtj/biolinks_vs_rtcga.xlsx?dl=0
  [5]: https://www.dropbox.com/s/d99djwqkuwxm1qn/RTCGA_v1.R?dl=0
  [6]: https://www.dropbox.com/s/lw2dkk9fb2wa4q8/TCGAbiolinks_v1.R?dl=0
Alright, I think I figured out what's going on.

The biggest difference in the results was caused by a missed point during pre-processing. Please see below:


`RTCGA` features an already preprocessed data frame with clinical parameters. In this data frame there is a `times` variable that combines both `days_to_death` and `days_to_last_follow_up`. That way you can provide `times` variable into the `survfit()` function and get the data in one step. In `TCGAbiolinks` data I needed to create this variable myself. Previously I overlooked this point, and I fitted the survival model only by using `days_to_death` variable. This effectively discards data from censored and alive patients skewing the results.


I also compared `legacy` and `harmonized` `miRNA` datasets by using TCGAbiolinks package. Results are shown below all side by side. As an example, two different parameters are plotted: `patient.gender` and `hsa-mir-155` expression. `harmonized` and `legacy` data uses the same clinical annotation file, that's why gender plots look the same. With this realization, `RTCGA` and `TCGAbiolinks` are giving comparable results. Whew! 


![KM comparison][1]


Thank you all for your insights!

Best,

Atakan
  [1]: https://s15.postimg.cc/91i5k280r/comparison.jpg