-
Notifications
You must be signed in to change notification settings - Fork 16
/
RiboseQC.R
141 lines (92 loc) · 7.24 KB
/
RiboseQC.R
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
## ----setup_vig, include=FALSE--------------------------------------------
knitr::opts_chunk$set(echo = TRUE)
## ----load_lib------------------------------------------------------------
library("RiboseQC")
tmp_dir <- tempdir()
file.path()
## ----create_dirs_hum-----------------------------------------------------
download.file("https://drive.google.com/uc?export=download&id=1n6UA5cSz6djx0dY_7sI57177T-kyZhhT",destfile = file.path(tmp_dir,"test_human.2bit"))
download.file("https://drive.google.com/uc?export=download&id=19Am9-iMEyB-AcIsVRdIqrF-BVdQYrcaI",destfile = file.path(tmp_dir,"test_human.gtf"))
## ----create_annot_hum----------------------------------------------------
prepare_annotation_files(annotation_directory = tmp_dir,
twobit_file = file.path(tmp_dir,"test_human.2bit"),
gtf_file = file.path(tmp_dir,"test_human.gtf"),scientific_name = "Human.test",
annotation_name = "genc25_22M",
export_bed_tables_TxDb = FALSE,
forge_BSgenome = TRUE,
create_TxDb = TRUE)
## ----load_hum------------------------------------------------------------
load_annotation("test_human.gtf_Rannot")
## ----gen_hum-------------------------------------------------------------
genome_seq[["chr22"]]
genome_seq[["chrM"]]
## ----gtf_hum_general_1---------------------------------------------------
GTF_annotation$exons_txs
## ----gtf_hum_general_2---------------------------------------------------
GTF_annotation$cds_txs
## ----gen_hum_cds---------------------------------------------------------
getSeq(genome_seq,GTF_annotation$cds_txs[[4]])
## ----gtf_hum_general_3---------------------------------------------------
GTF_annotation$start_stop_codons
## ----gtf_hum_general_4---------------------------------------------------
GTF_annotation$cds_txs_coords
## ----gtf_hum_general_5---------------------------------------------------
GTF_annotation$trann
## ----gtf_hum_general_6---------------------------------------------------
GTF_annotation$genetic_codes
getGeneticCode(GTF_annotation$genetic_codes["chr22","genetic_code"])
getGeneticCode(GTF_annotation$genetic_codes["chrM","genetic_code"])
## ----gtf_hum_general_7---------------------------------------------------
GTF_annotation$genome_package
## ----human_report_1------------------------------------------------------
if(!file.exists('test_human_hek.bam')) download.file("https://drive.google.com/uc?export=download&id=11PP5y2QH7si81rbEBJsOB-Lt3l_JowRW",destfile = file.path(tmp_dir,"test_human_hek.bam"))
## ----human_report_2------------------------------------------------------
RiboseQC_analysis(annotation_file=file.path(tmp_dir,"test_human.gtf_Rannot"),bam_files = file.path(tmp_dir,"test_human_hek.bam"),report_file = file.path(tmp_dir,"test_human_hek.html"),write_tmp_files = F)
## ----human_report_3,warning=FALSE, results='asis', fig.width=12, fig.height=10, dpi=120----
readRDS("test_human_hek.html_plots/rds/sample1_nucl_4_profiles_P_sites_metagene_subcodon_all")[[1]]
## ----yeast_1-------------------------------------------------------------
download.file("https://drive.google.com/uc?export=download&id=1ul_mqVimaavfhJs76R3TaOyTcUgwJvS2",destfile = file.path(tmp_dir,"test_yeast.2bit"))
download.file("https://drive.google.com/uc?export=download&id=1yffqKkfohGx9y7XRZO05vhjYUXgjiKBF",destfile = file.path(tmp_dir,"test_yeast.gtf"))
prepare_annotation_files(annotation_directory = tmp_dir,
twobit_file = file.path(tmp_dir,"test_yeast.2bit"),
gtf_file = file.path(tmp_dir,"test_yeast.gtf"),scientific_name = "yeast.test",
annotation_name = "yeast_custom",
export_bed_tables_TxDb = TRUE,
forge_BSgenome = TRUE,
create_TxDb = TRUE)
## ----yeast_2-------------------------------------------------------------
load_annotation("test_yeast.gtf_Rannot")
genome_seq[["II"]]
genome_seq[["Mito"]]
GTF_annotation$cds_txs
GTF_annotation$trann
## ----yeast_3-------------------------------------------------------------
download.file("https://drive.google.com/uc?export=download&id=1Nwlf7u-hHC_fX0fUJuJu_dr2WxwGazVJ",destfile = file.path(tmp_dir,"test_yeast_TCP_input.bam"))
download.file("https://drive.google.com/uc?export=download&id=18gQF2SxxsGQqVfhUw-wFo2dBt1mafiAB",destfile = file.path(tmp_dir,"test_yeast_TCP_RS.bam"))
download.file("https://drive.google.com/uc?export=download&id=1HBb58mOt1vmf58yqz-wkZzfxZhBl_tmO",destfile = file.path(tmp_dir,"test_yeast_TCP_SSU.bam"))
bam_filepath_y<-c("test_yeast_TCP_input.bam","test_yeast_TCP_RS.bam","test_yeast_TCP_SSU.bam")
RiboseQC_analysis(annotation_file=file.path(tmp_dir,"test_yeast.gtf_Rannot"),bam_files = bam_filepath_y,fast_mode = T,report_file = file.path(tmp_dir,"test_yeast_TCP.html"),sample_names = c("input","RS","SSU"),dest_names = file.path(tmp_dir,c("input","RS","SSU")),write_tmp_files = F)
## ----yeast_4, warning=FALSE, results='asis', fig.width=12, fig.height=10, dpi=120----
as_ggplot(readRDS("test_yeast_TCP.html_plots/rds/SSU_nucl_4_profiles_fivepr_metagene_subcodon_log2")[[1]])
as_ggplot(readRDS("test_yeast_TCP.html_plots/rds/RS_nucl_4_profiles_fivepr_metagene_subcodon_log2")[[1]])
## ----arab_1--------------------------------------------------------------
download.file("https://drive.google.com/uc?export=download&id=1gTxlR66fRX0GsHsWtJB7piLelM8jhZGu",destfile = file.path(tmp_dir,"test_arabidopsis.2bit"))
download.file("https://drive.google.com/uc?export=download&id=1jYbdvaGfVnMx38X1P85Ounyq9jGwRVgP",destfile = file.path(tmp_dir,"test_arabidopsis.gtf.gz"))
prepare_annotation_files(annotation_directory = tmp_dir,
twobit_file = file.path(tmp_dir,"test_arabidopsis.2bit"),
gtf_file = file.path(tmp_dir,"test_arabidopsis.gtf.gz"),scientific_name = "arabidopsis.test",
annotation_name = "araport11_custom",export_bed_tables_TxDb = TRUE,
forge_BSgenome = TRUE,
create_TxDb = TRUE)
## ----arab_2--------------------------------------------------------------
annotation="test_arabidopsis.gtf.gz_Rannot"
download.file("https://drive.google.com/uc?export=download&id=144PIK15iJrDSshMWsGbl9zQvT27Z1fzt",destfile = file.path(tmp_dir,"test_arabidopsis_root.bam"))
download.file("https://drive.google.com/uc?export=download&id=1VVdMi0ho0LiRCld6-3vzibtg_CyPSTj3",destfile = file.path(tmp_dir,"test_arabidopsis_shoot.bam"))
bam_filepath=file.path(tmp_dir,c("test_arabidopsis_root.bam","test_arabidopsis_shoot.bam"))
RiboseQC_analysis(annotation_file=file.path(annotation,bam_files = bam_filepath),fast_mode = TRUE,report_file = file.path(tmp_dir,"test_root_shoots.html"),dest_names = c("root","shoots"),sample_names = c("root","shoots"),write_tmp_files = FALSE)
## ----arab_3--------------------------------------------------------------
readRDS("test_root_shoots.html_plots/rds/all_samples_1_readlocdist")[[1]]
## ----arab_4, warning=FALSE, results='asis', fig.width=12, fig.height=16, dpi=150----
as_ggplot(readRDS("./test_root_shoots.html_plots/rds/shoots_ChrC_7_codonusage_positional_A-sites_per_codon_all_zscore")[[1]])
## ----end_sessinf---------------------------------------------------------
session_info()