---
title: "Annotating gene-type and gene-length for genes in reference-2020-A"
author: Anushka Gupta
date: March 29 2021
---

# This notebook has commands to annotate gene-length and gene-type for all the genes in reference-2020-A 
```{r setup}
knitr::opts_knit$set(root.dir = '~/Desktop/scRNA-seq_snRNA-seq_adipocyte_lineages/')
```

```{r}
# Using the same GTF file as used in cellranger while generating ref-2020A
gencode.v32.primary_assembly.annotation <- read.delim("additional_files/gencode.v32.primary_assembly.annotation.gtf", header=FALSE, comment.char="#")

# Subsetting the GTF to only have gene-level information, as well as gene start (V4) and end coordinate (V5), and gene-metadata (V9)
gencode=gencode.v32.primary_assembly.annotation[gencode.v32.primary_assembly.annotation$V3=="gene",c("V4","V5","V9")]
rm(gencode.v32.primary_assembly.annotation)

# Extracting ensemble gene ID, gene biotype and gene-length (V5-V4)
ensg.id=vector(mode = "character",length = nrow(gencode))
biotype=vector(mode = "character",length = nrow(gencode))
gene_l=vector(mode = "double",length = nrow(gencode))
for (i in 1:nrow(gencode)) { #change this
  x=strsplit(gencode$V9[i],split = " ")[[1]]
  ensg.id[i]=x[2]
  biotype[i]=x[4]
  gene_l[i]=gencode$V5[i]-gencode$V4[i]
}

# Reading in gene-level information provided by cellranger output
# This file can also be found in additional_files folder 
features.new <- read.delim("additional_files/ref2020a.tsv", header=FALSE, stringsAsFactors=FALSE)
features.new$biotype=NA
features.new$length=NA
for (i in 1:nrow(features.new)) {
  # Annotate features.tsv files with the gene-length and gene-type based on the 1st occurence of ENSG ID
  index=match(features.new$V1[i],ensg.id)
  features.new$biotype[i]= strsplit(biotype[index],split = ";")[[1]][1]
  features.new$length[i]=gene_l[index]
}
save(features.new,file = "rdata/features.ref2020A.RData")
```


