R Notes

查看当前位置

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}cn=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"
~~~
最后编辑于 ：
©著作权归作者所有,转载或内容合作请联系作者
  • 人面猴
    序言：七十年代末，一起剥皮案震惊了整个滨河市，随后出现的几起案子，更是在滨河造成了极大的恐慌，老刑警刘岩，带你破解...
    沈念sama阅读 212,080评论 6赞 493
  • 死咒
    序言：滨河连续发生了三起死亡事件，死亡现场离奇诡异，居然都是意外死亡，警方通过查阅死者的电脑和手机，发现死者居然都...
    沈念sama阅读 90,422评论 3赞 385
  • 救了他两次的神仙让他今天三更去死
    文/潘晓璐 我一进店门，熙熙楼的掌柜王于贵愁眉苦脸地迎上来，“玉大人，你说我怎么就摊上这事。” “怎么了？”我有些...
    开封第一讲书人阅读 157,630评论 0赞 348
  • 道士缉凶录：失踪的卖姜人
    文/不坏的土叔 我叫张陵，是天一观的道长。 经常有香客问我，道长，这世上最难降的妖魔是什么？ 我笑而不...
    开封第一讲书人阅读 56,554评论 1赞 284
  • 港岛之恋（遗憾婚礼）
    正文 为了忘掉前任，我火速办了婚礼，结果婚礼上，老公的妹妹穿的比我还像新娘。我一直安慰自己，他们只是感情好，可当我...
    茶点故事阅读 65,662评论 6赞 386
  • 恶毒庶女顶嫁案：这布局不是一般人想出来的
    文/花漫 我一把揭开白布。 她就那样静静地躺着，像睡着了一般。 火红的嫁衣衬着肌肤如雪。 梳的纹丝不乱的头发上，一...
    开封第一讲书人阅读 49,856评论 1赞 290
  • 城市分裂传说
    那天，我揣着相机与录音，去河边找鬼。 笑死，一个胖子当着我的面吹牛，可吹牛的内容都是我干的。 我是一名探鬼主播，决...
    沈念sama阅读 39,014评论 3赞 408
  • 双鸳鸯连环套：你想象不到人心有多黑
    文/苍兰香墨 我猛地睁开眼，长吁一口气：“原来是场噩梦啊……” “哼！你这毒妇竟也来了？” 一声冷哼从身侧响起，我...
    开封第一讲书人阅读 37,752评论 0赞 268
  • 万荣杀人案实录
    序言：老挝万荣一对情侣失踪，失踪者是张志新（化名）和其女友刘颖，没想到半个月后，有当地人在树林里发现了一具尸体，经...
    沈念sama阅读 44,212评论 1赞 303
  • 护林员之死
    正文 独居荒郊野岭守林人离奇死亡，尸身上长有42处带血的脓包…… 初始之章·张勋 以下内容为张勋视角 年9月15日...
    茶点故事阅读 36,541评论 2赞 327
  • 白月光启示录
    正文 我和宋清朗相恋三年，在试婚纱的时候发现自己被绿了。 大学时的朋友给我发了我未婚夫和他白月光在一起吃饭的照片。...
    茶点故事阅读 38,687评论 1赞 341
  • 活死人
    序言：一个原本活蹦乱跳的男人离奇死亡，死状恐怖，灵堂内的尸体忽然破棺而出，到底是诈尸还是另有隐情，我是刑警宁泽，带...
    沈念sama阅读 34,347评论 4赞 331
  • 日本核电站爆炸内幕
    正文 年R本政府宣布，位于F岛的核电站，受9级特大地震影响，放射性物质发生泄漏。R本人自食恶果不足惜，却给世界环境...
    茶点故事阅读 39,973评论 3赞 315
  • 男人毒药：我在死后第九天来索命
    文/蒙蒙 一、第九天 我趴在偏房一处隐蔽的房顶上张望。 院中可真热闹，春花似锦、人声如沸。这庄子的主人今日做“春日...
    开封第一讲书人阅读 30,777评论 0赞 21
  • 一桩弑父案，背后竟有这般阴谋
    文/苍兰香墨 我抬头看了看天上的太阳。三九已至，却和暖如春，着一层夹袄步出监牢的瞬间，已是汗流浃背。 一阵脚步声响...
    开封第一讲书人阅读 32,006评论 1赞 266
  • 情欲美人皮
    我被黑心中介骗来泰国打工， 没想到刚下飞机就差点儿被人妖公主榨干…… 1. 我叫王不留，地道东北人。 一个月前我还...
    沈念sama阅读 46,406评论 2赞 360
  • 代替公主和亲
    正文 我出身青楼，却偏偏与公主长得像，于是被迫代替她去往敌国和亲。 传闻我的和亲对象是个残疾皇子，可洞房花烛夜当晚...
    茶点故事阅读 43,576评论 2赞 349