R Notes

作者: 浩瀚之宇 | 来源:发表于2019-02-17 22:01 被阅读0次

查看当前位置

getwd()

getwd()
[1] "/media/albert/0E3711AB0E3711AB/Research/TCGA/GSE1323"

使用getwd函数来显示当前工作目录，使用setwd函数改变当前目录

setwd(/media/albert/0E3711AB0E3711AB/Research/TCGA/GSE1323/GSE1323")

列出当前文件夹里文件

list.files()
[1] "GSE1323" "gse1323.r"

包的安装源添加

#ChAMP包的安装
source("https://bioconductor.org/biocLite.R")
options(BioC_mirror="http://mirrors.ustc.edu.cn/bioc/")
biocLite("ChAMP")
如果是墙内，最好加上options选择中科大镜像，因为有好多依赖包需要安装

改变包的源

install.packages(c('europepmc', 'ggplotify', 'ggraph', 'ggridges'))
Installing packages into ‘/home/albert/R/x86_64-pc-linux-gnu-library/3.4’
(as ‘lib’ is unspecified)
Error in readRDS(dest) : error reading from connection
options(repos=structure(c(CRAN="https://mirrors.tuna.tsinghua.edu.cn/CRAN/")))

shell 环境下安装github上的包

Rscript -e 'devtools::install_local("enrichplot")'##已下载的包

R 环境安装github上的包

devtools::install_local("enrichplot")##已下载的包
devtools::install_github(c("GuangchuangYu/DOSE", "GuangchuangYu/clusterProfiler"))

流程化安装包

library(BiocManager)
deps <- c("pathview","clusterProfiler","ELMER", "DO.db","GO.db",
"ComplexHeatmap","EDASeq", "TCGAbiolinks","AnnotationHub",
"gaia","ChIPseeker","minet","BSgenome.Hsapiens.UCSC.hg19",
"MotifDb","MotIV", "rGADEM", "motifStack","RTCGAToolbox")
for(pkg in deps) if (!pkg %in% installed.packages()) install(pkg, dependencies = TRUE)
deps <- c("devtools","DT","pbapply","readr","circlize")
for(pkg in deps) if (!pkg %in% installed.packages()) install.packages(pkg,dependencies = TRUE)
devtools::install_github("BioinformaticsFMRP/TCGAWorkflowData")
devtools::install_github("BioinformaticsFMRP/TCGAWorkflow", dependencies = TRUE)

删除安装的某个包

remove.packages(GEOquery)

Load and Unload Packages in R

You load the fortunes package like this:

library(fortunes)

The detach() function will let you unload a package

detach(package:fortunes)

suppressMessage

不显示载入包提示信息：
suppressMessages(library(tidyverse))

缺失值

x[is.na(x)] = mean(x[!is.na(x)]) #填空值/补缺失值

矩阵 OR 向量

library(dplyr)  
dat <- read.csv("femaleMiceWeights.csv")  
control = filter(dat, Diet == "chow") %>% select(Bodyweight) %>% unlist #   filter(dat, Diet == "chow") %>% select(Bodyweight) 还是个矩阵  

> dat <- read.csv("femaleMiceWeights.csv")  
> dat  
   Diet Bodyweight
1  chow      21.51
2  chow      28.14
3  chow      24.04
4  chow      23.45
5  chow      23.68
6  chow      19.79
7  chow      28.40
8  chow      20.98
9  chow      22.51
10 chow      20.10
11 chow      26.91
12 chow      26.25
13   hf      25.71
14   hf      26.37
15   hf      22.80
16   hf      25.34
17   hf      24.97
18   hf      28.14
19   hf      29.58
20   hf      30.92
21   hf      34.02
22   hf      21.90
23   hf      31.53
24   hf      20.73

> dat[,2]
 [1] 21.51 28.14 24.04 23.45 23.68 19.79 28.40 20.98 22.51 20.10 26.91 26.25
[13] 25.71 26.37 22.80 25.34 24.97 28.14 29.58 30.92 34.02 21.90 31.53 20.73
> dat[,"Bodyweight"]
 [1] 21.51 28.14 24.04 23.45 23.68 19.79 28.40 20.98 22.51 20.10 26.91 26.25
[13] 25.71 26.37 22.80 25.34 24.97 28.14 29.58 30.92 34.02 21.90 31.53 20.73
> x=dat[,2]
> class(x)
[1] "numeric"
> control = filter(dat, Diet == "chow") %>% select(Bodyweight)
> control
   Bodyweight
1       21.51
2       28.14
3       24.04
4       23.45
5       23.68
6       19.79
7       28.40
8       20.98
9       22.51
10      20.10
11      26.91
12      26.25
> control[,1]
 [1] 21.51 28.14 24.04 23.45 23.68 19.79 28.40 20.98 22.51 20.10 26.91 26.25

dplyr

Extract the first, last or nth value from a vector

Description:

 These are straightforward wrappers around ‘[[’. The main advantage
 is that you can provide an optional secondary vector that defines
 the ordering, and provide a default value to use when the input is
 shorter than expected.

Usage:

 nth(x, n, order_by = NULL, default = default_missing(x))
 
 first(x, order_by = NULL, default = default_missing(x))
 
 last(x, order_by = NULL, default = default_missing(x))

Arguments:

   x: A vector
   n: For ‘nth_value()’, a single integer specifying the position.
      Negative integers index from the end (i.e. ‘-1L’ will return
      the last value in the vector).

      If a double is supplied, it will be silently truncated.

order_by: An optional vector used to determine the order

default: A default value to use if the position does not exist in the
input. This is guessed by default for base vectors, where a
missing value of the appropriate type is returned, and for
lists, where a ‘NULL’ is return.

      For more complicated objects, you'll need to supply this
      value. Make sure it is the same type as ‘x’.

Value:

 A single value. ‘[[’ is used to do the subsetting.

Examples:
x <- 1:10
y <- 10:1

 first(x)
 last(y)
 
 nth(x, 1)
 nth(x, 5)
 nth(x, -2)
 nth(x, 11)
 
 last(x)
 # Second argument provides optional ordering
 last(x, y)
 
 # These functions always return a single value
 first(integer())

NAs produced by integer overflow

dat[16,1]dat[16,2]
[1] NA
Warning message:
In dat[16, 1] * dat[16, 2] : NAs produced by integer overflow
as.numeric(dat[16,1]) as.numeric(dat[16,2])
[1] 16030361700
dat[16,1]
[1] 150
dat[16, 2]
[1] 106869078
150 * 106869078
[1] 16030361700

ggplot2

theme_classic

theme_classic(): 无栅格

library("ggplot2")
ggplot(data = mpg, aes(x = displ,y = hwy,color = class)) + geom_point( )

ggplot(data = mpg, aes(x = displ,y = hwy,color = class)) + geom_point( ) +  theme_classic()

显著性添加工具---ggsignif

ggsignif 是发表在github上的开源包，专门用于在box plot上添加显著性标签。
https://github.com/const-ae/ggsignif

安装稳定版本：
install.packages("ggsignif")
安装最新开发版本：
devtools::install_github("const-ae/ggsignif")

library(ggplot2)
library(ggsignif)    #载入ggsignif
#我们使用iris数据集作为演示，iris数据集Species作为分类标签，Species有3个类别（"versicolor"、"virginica"、"setosa"）
#Species的三组两两分别作差异性检验，提前设定好配对分析的list：
compaired_list <- list(c("versicolor", "virginica"), 
     c("versicolor","setosa"), 
     c("virginica","setosa"))
#绘制geom_boxplot()，选择好统计检验方法：
p <- ggplot(iris, aes(Species, Sepal.Width, fill = Species)) +
    geom_boxplot() +
    ylim(1.5, 6.5) +
    geom_signif(comparisons = compaired_list,
                step_increase = 0.3,
                map_signif_level = F,
                test = wilcox.test) +
    theme_classic()
p

##不显示p值，而显示***
p <- ggplot(iris, aes(Species, Sepal.Width, fill = Species)) +
    geom_boxplot() +
    ylim(1.5, 6.5) +
    geom_signif(comparisons = compaired_list,
                step_increase = 0.3,
                map_signif_level = T,
                test = wilcox.test) +
    theme_classic()
p

topTags: Table of the Top Differentially Expressed Tags

(edgeR)

Description

Extracts the top DE tags in a data frame for a given pair of groups, ranked by p-value or absolute log-fold change.

Usage

topTags(object, n=10L, adjust.method="BH", sort.by="PValue", p.value=1)

# generate raw counts from NB, create list object
y <- matrix(rnbinom(80,size=1,mu=10),nrow=20)
d <- DGEList(counts=y,group=rep(1:2,each=2),lib.size=rep(c(1000:1001),2))
rownames(d$counts) <- paste("gene",1:nrow(d$counts),sep=".")

# estimate common dispersion and find differences in expression
# here we demonstrate the 'exact' methods, but the use of topTags is
# the same for a GLM analysis
d <- estimateCommonDisp(d)
de <- exactTest(d)

# look at top 10
topTags(de)
# Can specify how many genes to view
tp <- topTags(de, n=15)
# Here we view top 15
tp
# Or ** order by fold change instead**
topTags(de,sort.by="logFC")

markdown 数理化公式

Codecogs 含数理化公式例子

Math and Statistics

Algebra

$\left(x-1\right)\left(x+3\right)$

$\sqrt{a^2+b^2}$

x = a_0 + \frac{1}{a_1 + \frac{1}{a_2 + \frac{1}{a_3 + a_4}}}

$x = a_0 + \frac{1}{a_1 + \frac{1}{a_2 + \frac{1}{a_3 + a_4}}}$

x = a_0 + \frac{1}{\displaystyle a_1 + \frac{1}{\displaystyle a_2 + \frac{1}{\displaystyle a_3 + a_4}}}

$x = a_0 + \frac{1}{\displaystyle a_1 + \frac{1}{\displaystyle a_2 + \frac{1}{\displaystyle a_3 + a_4}}}$

Trig

\cos^{-1}\theta
\sin^{-1}\theta

e^{i \theta}

\left(\frac{\pi}{2}-\theta \right )

Geometry

$\overrightarrow{AB}$

\overleftrightarrow{AB}

\widehat{AB}

$\Delta A B C$

Calculus

\frac{\partial y}{\partial x}

\frac{d}{dx}c^n=nx{n-1}

$\frac{d}{dx}e^{ax}=a\,e^{ax}$

\frac{d}{dx}\ln(x)=\frac{1}{x}

$\frac{d}{dx}\ln(x)=\frac{1}{x}$

\frac{d}{dx}\sin x=\cos x

$\sqrt[3]{\frac xy}$

$\sqrt[3]{\frac xy}$

$$ x = \dfrac{-b \pm \sqrt{b^2 - 4ac}}{2a} $$

$x = \dfrac{-b \pm \sqrt{b^2 - 4ac}}{2a}$

Matrices

$\begin{pmatrix} a_{11} & a_{12} & a_{13}\\ a_{21} & a_{22} & a_{23}\\ a_{31} & a_{32} & a_{33} \end{pmatrix}$

$\begin{pmatrix} a_{11} & \cdots & a_{1n}\\ \vdots & \ddots & \vdots\\ a_{m1} & \cdots & a_{mn} \end{pmatrix}$

sets

$\bigcup_{i=1}^{n}{X_i}$

$\bigcap_{i=1}^{n}{X_i}$

Statistics

${^n}C_r$

$\frac{n!}{r!(n-r)!}$

$\sum_{i=1}^{n}{X_i}$

$\sum_{i=1}^{n}{X_i^2}$

$X_1, \cdots,X_n$

$\frac{x-\mu}{\sigma}$

$\sum_{i=1}^{n}{(X_i - \overline{X})^2}$

Hypergeometric test: 超几何检验

$p=1-\sum_{k=0}^m \frac{\binom{K}{k}\binom{N-K}{n-k}}{\binom{N}{n}}$

$p=1-\sum_{k=0}^m \frac{\binom{K}{k}\binom{N-K}{n-k}}{\binom{N}{n}}$

-N: total number of genes (20000)
-K: total number of genes annotated with this term
-n: number of genes in the selected list
-k: number of genes in the list annotated with this term

$$t= \frac{\bar x_1 - \bar x_2}{\sqrt{\frac{s_1^2}{n1} + \frac{s_2^2}{n2}}}$$

$t= \frac{\bar x_1 - \bar x_2}{\sqrt{\frac{s_1^2}{n1} + \frac{s_2^2}{n2}}}$

Phisics

$\vec{F}=m\vec{a}$

$e=m c^2$

$\vec{F}=m \frac{d \vec{v}}{dt} + \vec{v}\frac{dm}{dt}$

$\oint \vec{F} \cdot d\vec{s}=0$

$\vec{F}_g=-F\frac{m_1 m_2}{r^2} \vec{e}_r$

Chemistry

$_{10}^{5}C^{16}$

$2H_2 + O_2 \xrightarrow{n,m}2H_2O$

$A\underset{b}{\overset{a}{\longleftrightarrow}}B$

$A\underset{0}{\overset{a}{\rightleftarrows}}B$

$A\underset{0^{\circ}C }{\overset{100^{\circ}C}{\rightleftarrows}}B$

Others

FPKM

FPKM是Fragments Per Kilobase per Million的缩写

$FPKM=\frac {10^6 \times n_f}{L \times N}$

$FPKM=\frac {10^6 \times n_f}{L \times N}$

其中， $n_f$ 是比对至目标基因的fragments的数量。 L是目标基因的外显子长度之和除以1000，单位是kb，不是bp；N是总有效比对至基因组的read数量。

表格

FDR(False Discovery Rate)

| |# not rejected <br>Not called |# rejected <br>Called |Total
------|----|----|----
# H 0 <br>Two groups similar |U |V |$m_0$
# H 1 <br>Two groups different |T |S |$m_1$

	# not rejected Not called	# rejected Called	Total
# H 0 Two groups similar	U	V	$m_0$
# H 1 Two groups different	T	S	$m_1$
Total	m-R	R	m

表格单元内字换行用<br>

R包的使用与实用技巧

TCGA数据ENSG转换为基因名(Symbol)

TCGA的RNAseq数据，是使用gencode进行基因注释的。因此，下载raw_counts后，基因名称是"ENSG"开头的ensemble号。ensemble号没有实际意义，而发表文章常用的是基因名称，经常阅读文献可以从基因名称猜测基因功能。因此，需要将ENSG转换为Symbol（基因名称）。

将ENSG转换Symbol，是常规操作，于是，首先选择找找已有的R包，而不是选择撸起柚子造轮子。


# 安装包

source("https://bioconductor.org/biocLite.R")

biocLite("AnnotationDbi")

biocLite("org.Hs.eg.db")

# 加载包

library("AnnotationDbi")

library("org.Hs.eg.db")

# for test

GENEID <- c(1,2,3,9,10)

ENSEMBL <- c("ENSG00000121410","ENSG00000175899","ENSG00000256069","ENSG00000171428","ENSG00000156006")

df <- data.frame(ENSEMBL,GENEID)

# ENSG转换Symbol

df$symbol <- mapIds(org.Hs.eg.db,

                     keys=ENSEMBL,

                     column="SYMBOL",

                     keytype="ENSEMBL",

                     multiVals="first")

#write.table "quote =F"和“quote =T(缺省)”的区别：
write.table(filt_set_annotated, file="filt_set_annotated.txt", row.names=F, sep="\t", quote =F)
~~~
ID  adj.P.Val   P.Value t   B   logFC   SPOT_ID SEQUENCE    Alias   source  chrom   strand  txStart txEnd   circRNA_type    best_transcript GeneSymbol
ASCRP002600 1.73691714077242e-11    5.00408280257107e-15    -35.5867037859303   24.2526719407331    -4.02735361814286   hsa_circRNA_102272  CCTAACCCAGGATGACGTCATGATGCCCCAGCGGGTGAGGCTGCAA  hsa_circ_0000059    circBase    chr1    +   40506473    40506637    intronic    NM_001105530    CAP1
~~~
write.table(filt_set_annotated, file="filt_set_annotated.txt", row.names=F, sep="\t", quote =T)

OR

write.table(filt_set_annotated, file="filt_set_annotated.txt", row.names=F, sep="\t")

~~~
"ID"    "adj.P.Val" "P.Value"   "t" "B" "logFC" "SPOT_ID"   "SEQUENCE"  "Alias" "source"    "chrom" "strand"    "txStart"   "txEnd" "circRNA_type"  "best_transcript"   "GeneSymbol"
"ASCRP002600"   1.73691714077242e-11    5.00408280257107e-15    -35.5867037859303   24.2526719407331    -4.02735361814286   "hsa_circRNA_102272"    "CCTAACCCAGGATGACGTCATGATGCCCCAGCGGGTGAGGCTGCAA"    "hsa_circ_0000059"  "circBase"  "chr1"  "+" "40506473"  "40506637"  "intronic"  "NM_001105530"  "CAP1"
~~~