R语言快速上手

参考
说明
- 安装程序包
数据类型
数学函数
判断
- if
- swith
- 与或非
循环
自定义函数
把source当模块用
绘图
- 函数
格式化输入输出
- cat输出信息到屏幕/文件
- print的输出会带很多信息，如上
文件读写
R调用Fortran
R调用Python
有趣
将R语言，变成Linux的小脚本使用
plot
ggplot2
- qplot
- ggplot
- 3d画图
R语言用于科学计算
- FFT
  - FFT
  - IFFT示例

R语言快速上手

参考

.Fortran（）返回任何结果
卡巴, 科弗, 高涛, 等. R 语言实战[J]. 2013.
dnorm, pnorm, qnorm与rnorm的区别
 R语言的流程控制
 R语言入门
 R 语言–接收命令行参数
 R里面有没有函数可以输出到已有文件末尾？
快速傅里叶变换 - MATLAB fft- MathWorks 中国
 【R语言数据导出txt】 write.table 函数用法

说明

运行脚本方法
1) R CMD [option] [xxx.R] [输出]
2) 与python和bash一样，在脚本里指明程序，添加可执行权限
```
#!/usr/bin/env Rscript
age <- c(1,2,3,4,5) 
weight <- c(2,3,4,5,6) 
plot(age,weight)
```
注释用#
同行书写，用;代表换行
变量名区分大小写,字母数字下划线
变量无需定义，随用随赋值，新创建，赋值为新数据=删除旧变量并重新定义赋值
R的数组和矩阵的元素类型不用相同
变量赋值，教程中是b <- 5，在R version 3.3.1可以直接b=5;b=c(1,2,3)
删除变量rm(变量名)
使用[” “]提取元素，有“，双引号
在脚本或程序中多次执行plot，每次的图会变成一个pdf的一页，程序执行结束后，有个总图的pdf保存在当前目录
使用source像模块一样调用
当前目录若存在.RData .Rhistory文件，好像会影响脚本运行结果

安装程序包

安装R语言

下载/编译安装清华-mirrors
可以下载安装包，也可以源码编译
homebrew安装
```
#MacOS也可以
brew install R
```

anaconda也可以安装

cndaqiang@girl:~/work/dfft-yimin/out$ conda install R

CRAN 进入R后安装，无需root

install.packages("inline")

指定清华镜像

options(repos=structure(c(CRAN="https://mirrors.tuna.tsinghua.edu.cn/CRAN/"))) 

Bioconductor

Bioconductor是大部分的生信相关的包存放的地方，网址https://www.bioconductor.org/ 安装使用代码如下

## try http:// if https:// URLs are not supported
source("https://bioconductor.org/biocLite.R")
biocLite("package.name")#将package.name替换为自己的要安装的包的名字即可。

GitHub

install.packages("devtools")
library(devtools)
install_github("user.name/package.name")#其中user.name是github作者的用户名。

查看程序包帮助

help(package = "sxtTools")
#如果想看某个具体的函数的信息，可以使用下面的代码：
?installBioc

数据类型

R的基本数据元素有数值，字符，逻辑型

基本数据元素

直接赋值，和数组就个c的区别

> a=5.0     #数值型
> a
[1] 5
> a="hello" #字符类型
> a
[1] "hello"
> a=T       #逻辑型
> a
[1] TRUE

数据类型转换

as.numeric(a) 转数值
as.character(b) 将数字b转化成字符格式
format(1.2345,digits = 3) 将数字转成字符格式，并设置有效位数为3
as.matrix(a)转矩阵
data.frame(Ave) 转数据框
as.numeric(as.vector(unlist(strsplit(Dipole,",")))) 组合操作，拆分字符串
options(digits=N),设置数据精度，屏幕显示的小数点位数

数组

数组元素可以是数值，字符，逻辑型混合在一起，不必相同类型

定义`c`

> a=4
> a
[1] 4
> a <- c(2,3)
> a
[1] 2 3

其他定义方式

> 1:4
[1] 1 2 3 4

调用`[ ]`

a[n]

矩阵

矩阵是由向量构成，因此元素可以是数值，字符，逻辑型混合在一起，不必相同类型

定义`matrix`

matrix(向量,[nrow=行数,ncol=列数][,byrow=T/F][,dimname=list(行名矩阵,列名矩阵)]  )

提取向量中的前nrow*ncol个元素填充矩阵
若不设置行数，默认为列向量
byrow=T向量按行填充，F按列填充，示例

> matrix(c(1,2,"hello",3),nrow=2,ncol=1)
     [,1]
[1,] "1"
[2,] "2"
> matrix(c(1,2,"hello",3),nrow=4,ncol=1)
     [,1]
[1,] "1"
[2,] "2"
[3,] "hello"
[4,] "3"
> matrix(1:4,nrow=2,ncol=2,byrow=F,dimname=list(c("hang1","hang2"),c("lie1","lie2")))
      lie1 lie2
hang1    1    3
hang2    2    4

调用`[ ]`

> b=matrix(1:4,nrow=2,ncol=2,byrow=F,dimname=list(c("hang1","hang2"),c("lie1","lie2")))
> b[1,1]
[1] 1
> b["hang1","lie2"]
[1] 3
> aa
     [,1] [,2] [,3] [,4]
[1,]    1    3    5    7
[2,]    2    4    6    8
> aa[1,2:3]
[1] 3 5

也支持像matlab一样，条件调用

> x=1:10
> x[x>5]
[1]  6  7  8  9 10

注意，a:b/2其实是c(a:b)/2, 应该写 a:(b/2),在提取矩阵时用, mydata[1:(n/2),2]

数组`array`

矩阵只能是二维的，数组与矩阵一样，区别是可以是高维的

定义`array`

示例

 b=array(1:8,c(2,2,2),dimname=list(c("hang1","hang2"),c("lie1","lie2"),c("matrix1","matrix2")))
 > b
, , matrix1

      lie1 lie2
hang1    1    3
hang2    2    4

, , matrix2

      lie1 lie2
hang1    5    7
hang2    6    8

特定数组

x=1:10 !等差数组
> seq(1,5,2)
[1] 1 3 5

提取`[ ]`

> b["hang2","lie1","matrix2"]
[1] 6

数据框

是将多个变量保存到一个数据框中
可以根据变量名从数据框中提取数据

定义`data.frame(col1,col2,col3...)`

每列元素要一样多

> col1=1:5
> col2=2:6
> col3=c("q","w","e","r","t")
> data.frame(col1,col2,col3)
  col1 col2 col3
1    1    2    q
2    2    3    w
3    3    4    e
4    4    5    r
5    5    6    t
#指定行名称
> data.frame(col1,col2,col3,row.names=col3)
  col1 col2 col3
q    1    2    q
w    2    3    w
e    3    4    e
r    4    5    r
t    5    6    t
#后期修改列名
dipoledate=data.frame(result$getDipolePython())
names(dipoledate)=c("Time","x","y","z")

调用`[ ]`

个人觉得$挺好用的

> b=data.frame(col1,col2,col3)
> b["col2"]
  col2
1    2
2    3
3    4
4    5
5    6
> b[1:3,2:3] #和matlab一样呀，[  [行],列]内只有一个数组则从列提取
  col2 col3
1    2    q
2    3    w
3    4    e
> b$col1
[1] 1 2 3 4 5
> summary(b$col1)
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max.
      1       2       3       3       4       5

也可以通过attach将数据框加到R的搜索路径中，就可以直接用

> rm(col1)
> col1
Error: object 'col1' not found
> attach(b)     #加入搜索路径
The following objects are masked _by_ .GlobalEnv:

    col2, col3

> col1
[1] 1 2 3 4 5
> detach(b)     #移除搜索路径
> col1
Error: object 'col1' not found

也可以使用with(数据框,{命令})

> with(b,{
+ col1
+ })
[1] 1 2 3 4 5

因子

因子可以理解为对已有数据的分类
一般含有字符型的数组，默认就是factor因子

> sex=factor(c("男","男","女","男"))
> sex
[1] 男 男 女 男
Levels: 女 男

指定设置分类顺序levels=c(...)

> sex=factor(c("男","男","女","男"),levels=c("男","女"))
> sex
[1] 男 男 女 男
Levels: 男 女

分类后，并不影响原数据,可以更好的用于数据分析

> sex[1]
[1] 男
Levels: 男 女

常将分类后的向量组成数据框，用于进一步分析

> sex
[1] 男     男     女     Unkown
Levels: Unkown 女 男
> num
[1] 12 14 12 15
Levels: 12 14 15
> name
[1] "Li" "Ya" "CQ" "LL"
> data=data.frame(name,sex,num)
> data
  name    sex num
1   Li     男  12
2   Ya     男  14
3   CQ     女  12
4   LL Unkown  15
> str(data)
'data.frame':   4 obs. of  3 variables:
 $ name: Factor w/ 4 levels "CQ","Li","LL",..: 2 4 1 3
 $ sex : Factor w/ 3 levels "Unkown","女",..: 3 3 2 1
 $ num : Factor w/ 3 levels "12","14","15": 1 2 1 3
> summary(data)
 name       sex    num
 CQ:1   Unkown:1   12:2
 Li:1   女    :1   14:1
 LL:1   男    :2   15:1
 Ya:1

列表

就像一个垃圾站，什么都能往里丢

定义

list(keywords1=...,keywords2=a,b,c...,keywords3....)

> list(one=num,two=name,sex)
$one
[1] 12 14 12 15
Levels: 12 14 15

$two
[1] "Li" "Ya" "CQ" "LL"

[[3]]
[1] 男     男     女     Unkown
Levels: Unkown 女 男

调用`listname[" keywords"]`

> lista["two"]
$two
[1] "Li" "Ya" "CQ" "LL"
> lista$two
[1] "Li" "Ya" "CQ" "LL"

其他

对变量注释`names`

> a=c("xiaoming",13)
> a
[1] "xiaoming" "13"
> names(a)[1]="Name"
> names(a)[2]="age"
> a
      Name        age
"xiaoming"       "13"

变量(数据对象)操作

class(变量名)查看变量类型

# 在拿到数据之前第一步要做的就是概览数据
head(df,n=10)#查看数据集的前十条记录
str(df)#查看数据集结构
tail(df)#查看数据集最后部分
sapply(df)#查看数据集类型及名称
summary(df)#查看数据集各变量描述统计

函数用途

length(object) 显示对象中元素/成分的数量
dim(object) 显示某个对象的维度
str(object) 显示某个对象的结构
class(object) 显示某个对象的类或类型
mode(object) 显示某个对象的模式
names(object) 显示某对象中各成分的名称
c(object, object,…) 将对象合并入一个向量 

cbind(object, object, …) 按列合并对象
rbind(object, object, …) 按行合并对象
Object 输出某个对象
head(object) 列出某个对象的开始部分
tail(object) 列出某个对象的最后部分
ls() 显示当前的对象列表
rm(object, object, …) 删除一个或更多个对象。语句rm(list = ls())
将删除当前工作环境中的几乎所有对象*
newobject <- edit(object) 编辑对象并另存为newobject
fix(object) 直接编辑对象，调用vi编辑

变量(数据对象)处理函数

字符

函数grep()、 sub()和strsplit() 搜索某个文本字符串（ fixed=TRUE）
搜索正则表达式（ fixed=FALSE，默认值为FALSE）
正则表达式可参考[gitbook文章重置] 正则表达式 grep sed awk

转义字符串

要使用特殊字符时，要在前面加\\，如

name=unlist(strsplit(as.character(mydata[id,1]),"\\|"))[3]

其他处理函数

矩阵操作

转置t(a)

> a=matrix(1:10,nrow=2)
> t(a)
     [,1] [,2]
[1,]    1    2
[2,]    3    4
[3,]    5    6
[4,]    7    8
[5,]    9   10

整合重构数据，暂略

数学函数

和matlab差不多啊

x=1:10
plot(x,sin(x))

也可以画x图

数学运算

abs, 取整, cos, sin, acos,asin,三角函数,ln=log(x,base=exp(1))

统计函数

mean,sum

概率函数

在R中，概率函数形如：

[dpqr][分布名称](参数表)

[dpqr]表示其所指分布的某一方面：

d = 密度函数（ density）!就是统计函数值
p = 分布函数（ distribution function） !统计函数的积分
q = 分位数函数（ quantile function）
r = 生成随机数（随机偏差）!生成符合分布的一组数x，x符合这一分布函数

以正态分布norm为例
其他分布名称

随机数种子

R每次计算伪随机数都会使用不同的随机数种子，我们可以指定随机数种子
这样可以重复出(自己/别人)的计算结果，如下
**注意：种子只生效一次，若仍使用该种子，要继续set.seed() **

> rnorm(10)
 [1]  1.2240818  0.3598138  0.4007715  0.1106827 -0.5558411  1.7869131
 [7]  0.4978505 -1.9666172  0.7013559 -0.4727914
> rnorm(10)
 [1] -1.0678237 -0.2179749 -1.0260044 -0.7288912 -0.6250393 -1.6866933
 [7]  0.8377870  0.1533731 -1.1381369  1.2538149
> set.seed(123)
> rnorm(10)
 [1] -0.56047565 -0.23017749  1.55870831  0.07050839  0.12928774  1.71506499
 [7]  0.46091621 -1.26506123 -0.68685285 -0.44566197
> set.seed(123)
> rnorm(10)
 [1] -0.56047565 -0.23017749  1.55870831  0.07050839  0.12928774  1.71506499
 [7]  0.46091621 -1.26506123 -0.68685285 -0.44566197

多元正态数据`mvrnorm` 略

将变量(矩阵，数组，数据框…)作为函数的输入参数

直接置入

> x=1:10
> sin(x)
 [1]  0.8414710  0.9092974  0.1411200 -0.7568025 -0.9589243 -0.2794155
 [7]  0.6569866  0.9893582  0.4121185 -0.5440211
> xx=matrix(1:10,nrow=2)
> sin(xx)
          [,1]       [,2]       [,3]      [,4]       [,5]
[1,] 0.8414710  0.1411200 -0.9589243 0.6569866  0.4121185
[2,] 0.9092974 -0.7568025 -0.2794155 0.9893582 -0.5440211
> mean(xx)
[1] 5.5

apply （数组，矩阵）某一维度

apply(变量名,维度，函数名)

> apply(xx,1,sin)
           [,1]       [,2]
[1,]  0.8414710  0.9092974
[2,]  0.1411200 -0.7568025
[3,] -0.9589243 -0.2794155
[4,]  0.6569866  0.9893582
[5,]  0.4121185 -0.5440211
> apply(xx,1,mean)
[1] 5 6

lapply,sapply (列表list)的某一维度，暂略

自编函数

判断

if

}后跟新的else

if (条件)
{ 
执行
}else if
{
执行    
}else{
执行 
}

示例

> if(a==10)
+ {
+ a=5
+ }
> if(a==5)
+ {
+ a=5
+ }else if(a==6)
+ {
+ a=6
+ }else
+ {a=10
+ }
> a
[1] 5

swith

swith(表达式,case1,case2,cse3)
#表达式为整数时，返回表达式对应的case值
#表达式为字符串时，case必须有名称，见下例，返回case值
#非整数非字符串，系统会转换成整数

好像就只能这么用

> a=2
> switch(a,"ming","hong","li")
[1] "hong"
> switch("a",a="ming",b="hong",c="li")
[1] "ming"

与或非

&、|、!、&&、||、xor

循环

for ( i in seq)

若省略{ }只循环for下面的一句，类似c

> for (i in 1:10)
+ {
+ y[i]=i
+ }

while (条件)

与for一样

> > while (i>5)
+ {
+ i=i-2
+ y=i
+ }

循环控制break和next

自定义函数

先定义后使用 示例

myfun=function(x,y=2,z=3)
{
a=x*y+z
return(a)
}
myfun(x=0)

调用其他文件中的函数，先

source("myfun.R")
y=myfun(x=0)

把source当模块用

source("myfun.R")

绘图

函数

plot 散点，点线
hist 直方图
boxplot
par 
...等等

线性，比例尺等也可以调控，都可以查表，或者单独做一个教程，此文不多讲，（估计没）有时间重开一个文章写绘图

pretty生成特定坐标，用于画图

pretty(x,n) 创建美观的分割点。通过选取n+1个等间距的取整值，将一个连续型变量x分割为n个区间
pretty(1:15)          # 0  2  4  6  8 10 12 14 16
pretty(1:15, h = 2)   # 0  5 10 15
pretty(1:15, n = 4)   # 0  5 10 15
pretty(1:15 * 2)      # 0  5 10 15 20 25 30
pretty(1:20)          # 0  5 10 15 20
pretty(1:20, n = 2)   # 0 10 20
pretty(1:20, n = 10)  # 0  2  4 ... 20

格式化输入输出

`cat`输出信息到屏幕/文件

#源码
#!/usr/bin/env Rscript
cat("the total number is:",3+5,"\n")
#sep指定各个间隔，file="输出到文件"
cat(1:4,"\n",sep="\t",file="catout") #默认清除源文件输入
cat(1:4,sep=c(",","!","@","\n"),file="catout", append = TRUE) # append = TRUE追加到末尾
#sep也可以指定每一个字符之间的间隔
#也可使用   sep=c(rep("\t",3),"\n") 
cat(1:4,sep=c(rep("\t",3),"\n"),file="catout", append = TRUE)
x=3+4
print(x)
#运行
cndaqiang@win10:~/code/R/Rcode> ./test.R
the total number is: 8
[1] 7
cndaqiang@win10:~/code/R/Rcode> cat catout
1       2       3       4
1,2!3@4
1       2       3       4

还可以这样

com=system("ls")
cat(com)
#结果
cndaqiang@win10:~/code/R/Rcode> ./test.R
catout    Cv.dat    myf.f90  myf.so           RCv.R       test      windowscode
cnq.R     ggplot.r  myf.mod  myfun.R          Rgraph.R    test.dat
cnq.Rout  input.R   myf.o    PErrorBar_CNQ.R  Rplots.pdf  test.R

`print`的输出会带很多信息，如上

文件读写

手输入`edit()`

new=edit()

然后填入c(1,2,3)即等价于new=c(1,2,3)
也可以在edit()括号中加入其他变量，进行修改,不影响原数据

> new=edit(c(2,3,4))

从文本文件读入

read.table("文件名"[,header=T/F,sep="间隔标识符",row.name="变量名"]) 
Cv=read.table("Cv.dat",skip=1) #skip跳过行号

header是否将数据文件的第一行识别为变量名,如

> read.table("my.dat")
 V1  V2  V3
1 num one two
2   1  11  21
3   2  12  22
4   3  13  23
> read.table("my.dat",header=T)
num one two
1   1  11  21
2   2  12  22
3   3  13  23

sep数据间隔，默认为一个或多个空格、制表符、换行符或回车符，默认就行

row.name=”变量名”，将变量名所在列作为行名

> read.table("my.dat",header=T,row.name="num")
one two
1  11  21
2  12  22
3  13  23

保存数据

write.csv(Dipole,"RDipole.csv",row.names =F)
#======保存图片
#
pdf("RDipole.pdf")
#png("filename.png")
#jpeg("filename.jpeg")
#bmp("filename.bmp")
#
#
#
#
for ( direct in 2:4)
{
plot(Dipole[,1],Dipole[,direct],type="l",main = "Dipole", xlab = names(Dipole)[1],  ylab = names(Dipole)[direct])
}
dev.off()

#控制画板
png(file="my.png",width=800,height=30000)

write.csv默认sep是`,`，其他具体参数同`write.table`

write.table (x, file ="", sep ="", row.names =TRUE, col.names =TRUE, quote =TRUE)

x：需要导出的数据
file：导出的文件路径
sep：分隔符，默认为空格（” “），也就是以空格为分割列
row.names：是否导出行序号，默认为TRUE，也就是导出行序号
col.names：是否导出列名，默认为TRUE，也就是导出列名
quote：字符串是否使用引号表示，默认为TRUE，也就是使用引号表示

貌似write.csv在输出data.frame时还是比不过write.table

> write.table(TDEFIELD_all,file="test",row.names =F, col.names =FALSE)
> write.csv(TDEFIELD_all,file="test",row.names =F, col.names =FALSE)
Warning message:
In write.csv(TDEFIELD_all, file = "test", row.names = F, col.names = FALSE) :
  不能修改'col.names'

其他

R调用Fortran

R会自动编译文件
myf.90文件

subroutine myfun(a,a2)
    REAL(8) :: a,a2 #Fortran必须使用双精度
    open(unit=20,file="test")
    write(20,*) a,a2
    close(20)
end subroutine myfun

subroutine ofun(a,b)
    REAL(8) :: a,b
    b=a**10
end subroutine ofun

cnq.R文件

#!/usr/bin/env Rscript
system("R CMD SHLIB myf.f90")  #编译动态库
#好像会自动根据时间戳，判断是否需要重新编译
dyn.load("myf.so") #导入动态库
#.Fortran("函数名",[索引1=]as.类型(输入数据),[索引2=]as.类型(输入数据).....)
#索引任意也可以没有
#计算结果为list型，使用[ ]或$提取计算结果
input=0
output=0
c=.Fortran("ofun",inp=as.double(2),outp=as.double(4))["outp"]
.Fortran("myfun",a=as.double(2),a2=as.double(c))

查看计算结果

cndaqiang@win10:/mnt/e/work/CODE/R/Rcode> cat test
   2.0000000000000000        1024.0000000000000

R调用Python

# 安装reticulate包
install.packages("reticulate")
# 加载reticulate包
library(reticulate)
> os=import("os")
> os$system("python --version")
Python 2.7.16 :: Anaconda, Inc.
[1] 0
result=import("pyramids.io.result")
dipo=result$getDipolePython()

 

有趣

`assign`将字符串转为变量

> x
 [1]  1  2  3  4  5  6  7  8  9 10
> assign('ww',x)#转为变量并赋值x
> ww
 [1]  1  2  3  4  5  6  7  8  9 10
 
#get()：返回与字符串同名的变量的值,并转成字符串
#assign():为字符串变量的字符串赋值
> get("xx")
 [1] "1" "1" "1" "1" "1" "1" "1" "1" "0" "0"
> 

`system("ls")`

调用系统命令system("系统命令")

com=system(c("cd ..;ls"))
#查看保存变量
cat(com)
#执行系统命令，并将结果保存到com
或者
> com=system(c("cd ..;ls"),intern=T)
> com
[1] "Ag55"        "Ag55.tar.gz" "graphenTDAP" "H2O"         "test-O3"

存储一些数据,函数,直接加载调用

> a='hello'
> myfun=function()
+ {
+ cat("hello R")
+ }
> save.image("cndaqiang.RData")
> q()

然后导入

# 导入数据
> load("cndaqiang.Rdata")
# 已经有了
> a
[1] "hello"
> myfun
function()
{
cat("hello R")
}
#查看变量空间有哪些内容
> ls.str()
a :  chr "hello"
myfun : function ()

将R语言，变成Linux的小脚本使用

在超算上编译安装R到个人目录

编译时候发现,R语言编译时有lapack，不错不错

wget https://mirrors.tuna.tsinghua.edu.cn/CRAN/src/base/R-3/R-3.6.0.tar.gz
tar xzvf R-3.6.0.tar.gz 
cd R-3.6.0/
module load  mpi/mvapich2/gnu/2.3b
./configure --prefix=/public/home/cndaqiang/soft/gcc-MVAPICH/R-3.6.0
make #不支持-j20
make install
PATH=/public/home/chendq/soft/gcc-MVAPICH/R-3.6.0/bin:$PATH
#--------------------------------------------
#Error ：configure: error: libcurl >= 7.22.0 library and headers are required with support for https
#参考https://www.chengzi520.com/?p=1912
wget https://curl.haxx.se/download/curl-7.57.0.tar.gz
tar zxvf curl-7.57.0.tar.gz
cd curl-7.57.0
./configure --prefix=/home/cndaqiang/software/gcc-4.8.5-openmpi-1.10.7/curl-7.57.0
make && make install 
# 在~/.bashrc添加：
#R
CURLDIR=/home/cndaqiang/software/gcc-4.8.5-openmpi-1.10.7/curl-7.57.0
export PATH=$CURLDIR/bin:$PATH
export LD_LIBRARY_PATH=$CURLDIR/lib:$LD_LIBRARY_PATH
#使环境变量生效：
source ~/.bashrc
#--------------------------------------------------
#WARNING：configure: WARNING: neither inconsolata.sty nor zi4.sty found: PDF vignettes and package manuals will not be rendered optimally 警告
#我们忽视了这个错误，有人解决如下
#https://my.oschina.net/u/2429108/blog/552977
3) 下载安装 inconsolata.sty
下载 ：wget http://mirrors.ctan.org/fonts/inconsolata.zip
解压 ：Unzip inconsolata.zip 
将文件拷贝到目录下：
        cp -Rfp inconsolata/* /usr/share/texmf/
 刷新sty ：
        mktexlsr

脚本

文件开头

#!/usr/bin/env Rscript
margs=commandArgs()
#参数
#print(margs)
#cndaqiang@win10:/mnt/e/work/CODE/R/Rcode> ./input.R 1 2 3
#[1] "/usr/lib64/R/bin/exec/R" "--slave"
#[3] "--no-restore"            "--file=./input.R"
#[5] "--args"                  "1"
#[7] "2"                       "3"
#linux中，第6个开始才是用户输入的参数
inpara=margs[6:length(margs)]
#常用包
#ggplot2
library(ggplot2)

快速画图

#!/usr/bin/env Rscript
##################################
# Author : cndaqiang             #
# Update : 2019-05-11            #
# Build  : 2019-05-11            #
# What   : 画x,y图               #
##################################
#头文件
#函数库
source("~/scripts/R-module/head.R")
#source(paste(CNQMODULEPATH,"xxx.R"))
datafile=inpara[1] #第一个输入参数为数据文件名
skipline=inpara[2] #第二个输入参数为是否跳过第一行
rplot=function(datafile,xlable="x",ylable="y",qtitle="",skipline=0)
{
#datafile="Cv.dat" #输入文件
#xlable="x" #x标签
#ylable="y" #y标签
#skipline=1 #跳过开头的几行
#读入数据
    mydata=read.table(datafile,skip=skipline)
    datacol=length(mydata[1,])
    datarow=length(mydata[,1])
    if (datacol == 1)
    {
        plot(mydata,type="b",main = qtitle, sub = "", xlab = xlable,  ylab = ylable)
    }else
    {

        ydata=mydata[,2:datacol] #
        xdata=mydata[,1] #第一列
        for ( i in 2:datacol )
        {
            y=mydata[,i]
            x=mydata[,1]
            plot(x,y,type="b",main = qtitle, sub = "", xlab = xlable,  ylab = ylable)
        }
    }
}

rplot(datafile,xlable=inpara[3],ylable=inpara[4],qtitle=inpara[5],skipline=skipline)

集成到makefile里,编译计算,画图,简直不要太爽

test:$(TEST)
	$(FC)  $(FFLAGS) -o test $(OBJS) $(LIBS)
	./test > result.dat
	./plot.r

调用shell的`grep`组合操作

提取矩阵

Dipole=paste("grep 'Electric dipole (a.u.)' ",result,"|grep 'TDAP' |awk '{printf \"%f,%f,%f,\",  $(NF-2),$(NF-1),$NF }' ")
Dipole=system(Dipole,intern=T)
Dipole=as.numeric(as.vector(unlist(strsplit(Dipole,",")))) #分离后，默认是list，转为向量和数值型
Dipole=c(Dipole)
Dipole=matrix(Dipole,ncol=3,byrow=T)

读入fdf格式

#======读入时间
MD.LengthTimeStep=paste("grep -i '^MD.LengthTimeStep'",result,"|awk '{ print $2 }'")
#此处默认单位是fs了，不对fs进行判断
MD.LengthTimeStep=as.numeric(system(MD.LengthTimeStep,intern=T))
MD.FinalTimeStep=paste("grep -i '^MD.FinalTimeStep'",result,"|awk '{ print $2 }'")
MD.FinalTimeStep=as.numeric(system(MD.FinalTimeStep,intern=T))

读入fdf block

#行号
phyHang=paste("grep -ni '^%block myblock'",input.fdf,"|awk -F: '{ print $1 }'")
phyHang=paste("hang=$(",phyHang,")")
phyHang=paste(phyHang,"hang=$(echo ${hang}+1|bc)",sep=";")
#读入
phy=paste(phyHang,";sed -n ${hang}p",input.fdf)
phy=system(phy,intern=T)
#
phy=as.vector(unlist(strsplit(phy," +"))) #以空格分开，strsplit生成的是list，要转为向量
phy=c(as.numeric(phy))

phydata=phy[1]       #有物理含义
phydata=phy[2]   #有物理含义
phydata=phy[3]       #有物理含义
phydata=phy[4]           #有物理含义
phydata=phy[5]           #有物理含义

plot

使用plot进行画图时，数据框的一列，一定要写成mydata[,lie]，把一列提成数组，那个,不能少示例，

pdf(file="RCurrentFFT.pdf")
for ( direct in 2:colnum)
{
plot(myfft[1:myn/2,1],myfft[1:myn/2,direct],type="l",main = "Current FFT", xlab = paste(names(myfft)[1],"FFT"), ylab = paste(names(myfft)[direct],"FFT") )
#----下面为加text到图中
minfft=max(myfft[1:myn/2,direct])/10
peak=(myfft[2:(myn/2),direct] > myfft[1:(myn/2-1),direct]) & (myfft[2:(myn/2),direct] > myfft[3:(myn/2+1),direct]) & ( myfft[2:(myn/2),direct] > minfft ) 
peak=c(T,peak,(1:(myn-2) < 0))
text(myfft[peak,1],myfft[peak,direct],format(myfft[peak,1],digits = 3),cex=1,pos=1,col="red")
#x,y,显示的内容
}
dev.off()

ggplot2

代码的发展，很多函数都变了，注意，有的教程给的函数不能用

qplot

qplot(数据框索引,数据框索引,data=数据框)
> price=1:10
> num=sin(price)
> mydata=data.frame(price,num)
qplot(price,num,data=mydata)
#也可以使用数学函数处理
qplot(price,log(x),data=mydata)
#也可以是数据框内多组数据的计算结果
qplot(price,log(x)*num,data=mydata)

根据分类，确定颜色和尺寸

> qplot(price,num,data=data,colour=good,shape=xx,size=num,alpha=num/10)
#shape=字符型，不然报错
#color字符数子都行，图形颜色
#size指定大小
#alpha透明度[0-1]之间
> mydata
   price        num xx good
1      1  0.8414710  1    0
2      2  0.9092974  1    0
3      3  0.1411200  1    0
4      4 -0.7568025  1    0
5      5 -0.9589243  1    0
6      6 -0.2794155  1    1
7      7  0.6569866  1    1
8      8  0.9893582  1    1
9      9  0.4121185  0    1
10    10 -0.5440211  0    1

geom指定形状

qplot(price,num,data=data,geom="line")
#geom="path","line"连线
#geom="point" 默认散点
#"smooth"平滑，使用span,method等继续控制，
#qplot(num,price,data=data,geom=c("smooth"),span=0.5,method="gam")
#也可以组合多种图形
qplot(price,num,data=data,geom=c("point","line"))

保存

ggsave("my.png")

貌似只保存最近的图

ggplot

#画图
png(outfile)
mypicture=ggplot(mydata,aes(x=xdata,y=ydata,color=orderdata))+geom_point()
#加水平垂直参考线(阈值线),由输入参数Qylimit和Qxlimit控制
mypicture=mypicture + geom_hline(yintercept=Qylimit,linetype=4)
mypicture=mypicture + geom_vline(xintercept=c(-1*Qxlimit,Qxlimit),linetype=4)
#设置x，y上下限
mypicture=mypicture + xlim(Lxlimit,Hxlimit)
mypicture=mypicture + ylim(Lylimit,Hylimit)
#设置点的颜色
mypicture=mypicture + scale_color_manual(values =c("Up" = "red", "Down" = "blue", "Normal" = "grey"))
#x，y轴label
#expression是为了输出上下标
mypicture=mypicture + labs(x=expression(log[2](FC)),y=expression(-log[10](FDR)),color="significant")

mypicture   
dev.off()

3d画图

zypper install Mesa-libGL-devel

R语言用于科学计算

FFT

数学课本f(t)-c(w)
代码t=t0+n*dt, f&n->FFT->C(n)
自己将n转为f(n),T(n)

#!/usr/bin/env Rscript
##################################
# Author : cndaqiang             #
# Update : 2019-05-20            #
# Build  : 2019-05-20            #
# What   : FFT分析               #
##################################

fn1=5;fn2=10;fn3=100;fn4=233 #频率
dt=0.002 #采样间隔
t=seq(0,1,dt) #样本点
n=length(t)  #采样数
y=sin(2*pi*fn1*t)+3*cos(2*pi*fn2*t)+4*sin(2*pi*fn3*t)+5*sin(2*pi*fn4*t) #采样值
#plot(x,y)
z=abs(fft(y)/n)  #频谱 
f=(1:n-1)/n/dt  #频率
png("fft.png")
plot(f[1:n/2],z[1:n/2],type="l")
dev.off()

结果图
图中的峰值比例不对，加密N点减少dt后就准了，就不补图了
由于输入是实数,这里f=(1:n-1)/n/dt #频率应该变为f=(1:n-1)/(n/2)/dt #频率

IFFT示例

FFT的结果需要除以N才能得到正确的强度

> a=1000
> x=1:a;x=x*2*pi/a;y=exp(1i*x);z=fft(y);  z1=fft(fft(y),inverse = TRUE);  z2=fft(fft(y)/a,inverse = TRUE)
> max(abs(z))
[1] 1000
> max(abs(z1))
[1] 1000
> max(abs(zz))
[1] 10000
> max(abs(z2))
[1] 1
> 

本文首发于我的博客@cndaqiang.
本博客所有文章除特别声明外，均采用 CC BY-SA 4.0 协议，转载请注明出处！

R语言快速上手

参考

说明

安装程序包

安装R语言

CRAN 进入R后安装，无需root

Bioconductor

GitHub

查看程序包帮助

数据类型

基本数据元素

数据类型转换

数组

定义c

调用[ ]

矩阵

定义matrix

调用[ ]

数组array

定义array

特定数组

提取[ ]

数据框

定义data.frame(col1,col2,col3...)

调用[ ]

因子

列表

定义

调用listname[" keywords"]

其他

对变量注释names

变量(数据对象)操作

变量(数据对象)处理函数

字符

转义字符串

其他处理函数

矩阵操作

整合重构数据，暂略

数学函数

数学运算

统计函数

概率函数

随机数种子

多元正态数据mvrnorm 略

将变量(矩阵，数组，数据框…)作为函数的输入参数

直接置入

apply （数组，矩阵）某一维度

lapply,sapply (列表list)的某一维度，暂略

自编函数

判断

if

swith

与或非

循环

for ( i in seq)

while (条件)

循环控制break和next

自定义函数

把source当模块用

绘图

函数

格式化输入输出

cat输出信息到屏幕/文件

print的输出会带很多信息，如上

文件读写

手输入edit()

从文本文件读入

保存数据

write.csv默认sep是,，其他具体参数同write.table

其他

R调用Fortran

R调用Python

有趣

assign将字符串转为变量

system("ls")

存储一些数据,函数,直接加载调用

将R语言，变成Linux的小脚本使用

在超算上编译安装R到个人目录

脚本

集成到makefile里,编译计算,画图,简直不要太爽

定义`c`

调用`[ ]`

定义`matrix`

调用`[ ]`

数组`array`

定义`array`

提取`[ ]`

定义`data.frame(col1,col2,col3...)`

调用`[ ]`

调用`listname[" keywords"]`

对变量注释`names`

多元正态数据`mvrnorm` 略

`cat`输出信息到屏幕/文件

`print`的输出会带很多信息，如上

手输入`edit()`

write.csv默认sep是`,`，其他具体参数同`write.table`

`assign`将字符串转为变量

`system("ls")`

调用shell的`grep`组合操作