SeqQuery/seq_query_data_processing.R at master · ginnyintifa/SeqQuery · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
# in silico translation

## data processing script
library(dplyr)
library(magrittr)
library(data.table)

tgf = readLines("/data/ginny/GENCODE/gencode.v27.basic.annotation.gtf")

### only extract the lines with "protein_coding"


tgf_pc = grep("protein_coding", tgf, value = T)


### make the file tab delimited


tgf_df = rbindlist(lapply(1:length(tgf_pc), function(x) {

  if(x%%100000 == 0)
    cat(x,"\n")

  ss = strsplit(tgf_pc[x], split = "\t", fixed = T)

  dss = data.frame(t(unlist(ss)), stringsAsFactors = F)

  return(dss)

}))


colnames(tgf_df) = c("chr","source","type", "start_position","end_position",
                     "score","strand", "frame", "description")


### only get CDS

tgf_df_CDS = tgf_df %>%
  dplyr::filter(type == "CDS")


parsed_desc = rbindlist(lapply(1:nrow(tgf_df_CDS),function(x) {

  this_desc = tgf_df_CDS$description[x]

  getit = unlist(strsplit(this_desc, split = ";", fixed = T))

  gene_id  = unlist(strsplit(getit[1], split = "\"",fixed = T))[2]
  transcript_id = unlist(strsplit(getit[2], split = "\"",fixed = T))[2]
  gene_type  = unlist(strsplit(getit[3], split = "\"",fixed = T))[2]
  gene_name = unlist(strsplit(getit[4], split = "\"",fixed = T))[2]
  transcript_type = unlist(strsplit(getit[5], split = "\"",fixed = T))[2]
  transcript_name = unlist(strsplit(getit[6], split = "\"",fixed = T))[2]
  exon_number = as.integer(unlist(strsplit(getit[7], split = " ",fixed = T))[3])
  exon_id  = unlist(strsplit(getit[8], split = "\"",fixed = T))[2]
  level = as.integer(unlist(strsplit(getit[9], split = " ",fixed = T))[3])
  protein_id = unlist(strsplit(getit[10], split = "\"",fixed = T))[2]
  havana_gene = unlist(strsplit(getit[14], split = "\"",fixed = T))[2]
  havana_transcript = unlist(strsplit(getit[15], split = "\"",fixed = T))[2]

  get_df = data.frame(gene_id, transcript_id, gene_type, gene_name,
                      transcript_type, transcript_name, exon_number, exon_id,
                      level, protein_id, havana_gene, havana_transcript, stringsAsFactors = F)

  if(x %% 10000 == 0)
    cat(x, "\n")

  return(get_df)
} ))


tgf_df_CDS$start_position = as.integer(tgf_df_CDS$start_position)
tgf_df_CDS$end_position = as.integer(tgf_df_CDS$end_position)


tgf_df_CDS_desc = cbind(tgf_df_CDS, parsed_desc)

tgf_df_CDS_desc = tgf_df_CDS_desc[,-9]

saveRDS(tgf_df_CDS_desc,file = "tgf_df_CDS.Rds")


aa = c("F","F","L","L",
       "S","S","S","S",
       "Y","Y","STOP","STOP",
       "C","C","STOP","W",
       "L","L","L","L",
       "P","P","P","P",
       "H","H","Q","Q",
       "R","R","R","R",
       "I","I","I","M",
       "T","T","T","T",
       "N","N","K","K",
       "S","S","R","R",
       "V","V","V","V",
       "A","A","A","A",
       "D","D","E","E",
       "G","G","G","G")


triple = c("UUU","UUC", "UUA", "UUG",
           "UCU","UCC","UCA","UCG",
           "UAU","UAC","UAA","UAG",
           "UGU","UGC","UGA","UGG",
           "CUU","CUC","CUA","CUG",
           "CCU","CCC","CCA","CCG",
           "CAU","CAC","CAA","CAG",
           "CGU","CGC","CGA","CGG",
           "AUU","AUC","AUA","AUG",
           "ACU","ACC","ACA","ACG",
           "AAU","AAC","AAA","AAG",
           "AGU","AGC","AGA","AGG",
           "GUU","GUC","GUA","GUG",
           "GCU","GCC","GCA","GCG",
           "GAU","GAC","GAA","GAG",
           "GGU","GGC","GGA","GGG")

triplet_dic = structure(aa, names = triple)


saveRDS(triplet_dic, file = "codon_dictionary.Rds")


wf = readLines("all_file/GRCh38.p10.genome.fa")

sig = grepl(">", wf)
starts = which(sig==T)


for(i in 1:24)
{
  startline = starts[i]+1
  endline = starts[i+1]-1
  seq = paste0(wf[startline:endline], collapse = "")
  cat(i, "\n")

  saveRDS(seq, file = paste0("chr", i,"_seq.Rds"))

}