「不定期更新」常用但又记不住的生信分析小技巧(code)

Section 1 R

  • 批量导入Table 储存为数据框: 比如有大量的数据,逗号分隔符的例如.csv,制表符分隔符的例如 .txt 部分 .xls,一个一个read太慢,需要批量导入
setwd("DirofYourData")
file_names <- list.files(pattern = ".xls")# 将路径下.xls的文件名储存到file_names这个变量中
for (i in 1:length(file_names)) { # for循环开始length(file_names)指的是该向量的长度,可以理解为从第一个到最后一个
  name<-gsub(".xls","",file_names[i])# file_names是一个向量用[]来索引那么file_names[1]表示上面赋值后第一个.xls文件名类似"lncRNA.fiber.xls",当i = 1时, 就变成gsub(".xls","",lncRNA.fiber.xls), 那么就是将.xls前面的文件名截取出来赋值给name。
  assign(name,read.table(file_names[i],sep = "\t",header = TRUE,quote = "",stringsAsFactors = FALSE))# assign是循环时给变量赋值,那么这段命令的意思就是依次从第一个到最后一个读入数据并赋值给name。
}

Reference:
Source: 百度, by 爸爸妈妈讲故事
Source: 简书,by 郑宝童

  • 批量导入数据并保存到list中
    Dec 10,2019 更新
    上面的方法如果数据很多,会有点不方便,下面的方法可以批量导入list中:
# 首先设置一下目标文件的路径
path.raw <- "~/heterosis_project/expression_model/03.enrichment/GO Enrichment/"
# 这里也是用pattern定义下你要读入的文件信息(文件名中共有的标志)
fileNames.raw <- dir(path.raw, pattern = "*.xls") 
filePath.raw <- sapply(fileNames.raw, function(x){ 
  paste(path.raw,x,sep='/')})   
data.raw <- lapply(filePath.raw, function(x){
  read.table(x, sep = "\t",header = TRUE,quote = "",stringsAsFactors = FALSE)}) 

Reference:
😹😹,暂时忘了,以后补上
最近改良了一下写成了function,我载考虑把我常用的function打成一个R package,不过目前没时间折腾,先source用
代码如下:
我习惯把数据都整理成Tab分割的txt,如果有人习惯csv或者其他,可以对应修改,也可以写个判断,加上参数type。

BatchReadTable = function(path, pattern,sep = "\t", header = TRUE, quote = "", stringsAsFactors = FALSE){
  fileNames.raw <- dir(path, pattern = pattern) 
  filePath.raw <- sapply(fileNames.raw, function(x){ 
    paste(path,x,sep='/')})
  data.raw <- lapply(filePath.raw, function(x){
    read.table(x, sep = sep,header = header,quote = "",stringsAsFactors = FALSE)}) 
  x = list(Name = fileNames.raw,
           file = data.raw)
  return(x)
}
#' Usage: 
path = "/filepath/"
pattern = "*.xls"
filelist = BatchReadTable(path = path,pattern = pattern) 
name = filelist$name; file = filelist$file
## name就是文件名,后续可以用gsub简化 例如name = filelist$name %>% gsub("*.xls","",.)
## file就是和上面那个代码一样输出的data.raw
  • 批量删除环境变量中的数据:做完数据清洗往往一些很大的原始数据还占据在环境变量中,所以想要删除,清空的命令rm(list = ls()),但是需要删除一些特定的含有共同字符的怎么办?刚开始想到通配符* rm(lncRNA.*),当然这个肯定报错,然后在高人指点下完成了,在此感谢BMK原老师🌹
rm(list=ls(pattern="lncRNA."))# 比如我之前导入的rawdata是lncRNA的FPKM,批量删除那么就让pattern="lncRNA"

An other solutionrm(list=ls(pattern="temp")), remove all objects matching the pattern.

这里的解释只要你的变量名中包含"lncRNA"都会被删除,所以这个是变通的,必须细心观察不然会遭到误伤😂
Reference:
Source: stack overflow

  • select()提取数据框感兴趣的列
    数据分析时通常需要从一个大的dataframe中提取感兴趣的列,或者数据清洗时需要删掉一些没有用的列,通过select()函数直接提取
# 比如我想把fpkm表达矩阵中叶片的fpkm值提取出来,那么
leaf.lncRNA = select(lncRNA.leaf, contains("leaf"))# 包含leaf的
# 再比如我再linux下paste差异列表的时候有些材料可能多次出现,但是导入R之后对于重复的colname,R会自动的在后面加.1.2.3...所以想要去掉这些重复的只取唯一的:
leaf.lncRNA = select(lncRNA.leaf, ends_with("FPKM"))以FPKM结尾的

Reference
墙裂推荐 by 简书 王诗翔

  • left_join做近似vlookup匹配
    left_join从字面上看就是向左融合,要求source中药匹配的列是唯一的,应该不支持一对多
result <- left_join(query, source, by = "comm"
# 这里要求两个dataframe中要匹配,比如要找gene_id对应的信息必须保证两个dataframe的gene_id那一列的colname都是gene_id.如果不一样修改也比较简单
query$id = query$x
source$id = source$y
# x和y分别代表query和source中gene id 那列的colname

其实对于两个dataframe的join操作远不止left_join只是我用的次数多,更多诸如 full_join, inner_join, left_join和right_join参考
Reference:
Source: 简书 WortJohn

  • counts TPM 和 FPKM互相转换
    这个好多人估计都有需求,直接搬运的,这些function建议自己建一个R project,然后存到环境变量中,可以直接用用法也是直接apply(inputmatrix,2,fpkmToTpm)
countToTpm <- function(counts, effLen)
{
    rate <- log(counts) - log(effLen)
    denom <- log(sum(exp(rate)))
    exp(rate - denom + log(1e6))
}
 
countToFpkm <- function(counts, effLen)
{
    N <- sum(counts)
    exp( log(counts) + log(1e9) - log(effLen) - log(N) )
}
 
fpkmToTpm <- function(fpkm)
{
    exp(log(fpkm) - log(sum(fpkm)) + log(1e6))
}
 
countToEffCounts <- function(counts, len, effLen)
{
    counts * (len / effLen)
}

Reference:
What the FPKM?
Source:简书 by 刘小泽

  • GOBG的格式转换
    在做GO富集分析时,许多工具需要提供 customized GO background, 本质上就是一个Gene_ID等代表这个基因或者转录本的符号对应GOterm,当然如果了解Gene ontology的定义及原理的话,很明显会知道在全基因层面(background)这种对应关系一定是多对多的;所以在格式上也出现了很多形式,第一种是gene_id是unique,每个gene_id只出现一次,但是在对应的GO_id部分,如果对应了多个GO_id,与之对应的GO_id会以空格或者逗号,分号等多种分割形式并排排列,如果在excel中可以表述为这些GO_id在同一单元格;此外还有一种形式是一个geneID对应一个GO_id,这时的Gene_id是重复的:例如

Type1

GeneID GO_ID
Gene1 GO1,GO2,GO3
Gene2 GO2
Gene3 GO10,GO11,GO12,GO100

Tpye2

GeneID GO_ID
Gene1 GO1
Gene1 GO2
Gene1 GO3
Gene2 GO2
Gene3 GO10
Gene3 GO11
Gene3 GO12
Gene3 GO100
###############################
#  convert type1 --> type2
###############################
# 读取第一种类型的文件
x = read.table("~/Dropbox/tools/inputquery2go.txt", 
               header = F, sep = "\t") 
head(x)
# install.packages("tidyverse")
library(splitstackshape)
# 由于我的文件中goid之间是空格分割的,所以这里在选择分割符的时候使用空格
y = cSplit(x, "V2", " ", "long")
# 如果成功列数肯定会成倍增加可以用dim查看,如果Rstudio直接在环境变量窗口检查,也可以用str(y)检查
dim(y)
###############################
#  convert type2 --> type1
###############################
# 原理暂且不清楚
z = aggregate(y, by = list(y$V1), c) %>% transmute(Gene = .$Group.1, GOid = .$V2)
# %>%的意思就是把左键的值发送给右键的表达式,并作为右键表达式函数的第一个参数。 类似shell的管道符 |
dim(x)
[1] 41939     2
dim(y)
[1] 119867      2
dim(z)
[1] 41939     2

其实我们可以发现上面用R做的type2-->type1的结果文件中其实每行都是一个list,在输出的时候比较麻烦,所以万能的excel又来了,这是你可能需要一个多核心,大内存的电脑,用excel的function实现了:

# 前两列就如同上面type2的格式排列就好
# 第三列(C)
=IF(COUNTIF(A$1:A1,A1)>1,"",A1)
# 第四列(D)
=IF(COUNTIF(A$1:A1,A1)>1,"",SUBSTITUTE(SUBSTITUTE(PHONETIC(OFFSET(A1,,,COUNTIF(A:A,A1),2)),A1,","),",","",1))

Reference:
Source: 简书, by 二货潜 😂每次引用这位老兄的文章感觉像是在骂人,强烈推荐,都是干货。👍👍
Source:百度知道,by 答主:jjchangyuan

Oct 22, 2019补充

  • RMarkdown中par函数layout会报错
    Rmarkdown是个不错的工具,可以直接将跑的结果以报告的形式整理出来,可以采用markdown的语法进行排版,一些HTML的语法也可以用,比如用<font siez=3 color = red></font>等来调整字体颜色等,但是发现一件蛋疼的事,把Rscript下的代码直接粘贴过来会发现有些会报错,比如用par函数调整layout,会报各种错误,后来一顿google发现需要在画图的时候用{}将画图的代码包起来,就不会报错了。
par(mfrow=c(1,2))
col.group <- group
levels(col.group) <-  brewer.pal(nlevels(col.group), "Set1")
col.group <- as.character(col.group)
col.type <- type
levels(col.type) <-  brewer.pal(nlevels(col.type), "Set2")
col.type <- as.character(col.type)
{plotMDS(lcpm, labels=group, col=col.group)
title(main="A. Sample groups")
plotMDS(lcpm, labels=type, col=col.type, dim=c(3,4))
title(main="B. Sample types")} #最好用{}把这里画图的代码包起来。
# glimma包看动态图
glMDSPlot(lcpm, labels=paste(group, type, sep="_"), 
          groups=y$samples[,c(1,4)], launch=T)
y$samples$lib.size

Dec 20, 2019补充

  • R删除向量中开头或者结尾的几个字符:
    做转录组的时候用的GID是isoform ID,造成在与其他databaseID不统一,存在.1,.2isoform的标记,其实这里需要做的就是删除结尾的两个字符就可以,这里用sub函数:
Gene ID Transcript ID
Gh_A01G000100 Gh_A01G000100.1
## 这里正则表达式..$表示结尾处两个字符这样去掉了字符结尾的后两位
GID = sub(pattern = "..$",replacement = "",x = GID2SwAc$V1)
## 同理如果要删掉开头的几个字符需要用^..
SwAc = sub(pattern = "^..",replacement = "",x = GID2SwAc$V2)

Section 2 Linux

  • Blast本地数据库构建
    这玩意不经常用老是记不住,每次构建都得help😂,老年人记性不好,下次直接翻简书
makeblastdb -in input_file -parse_seqids -hash_index -dbtype nucl
# makeblastdb: 就是make blast database
# -in 要建立本地blast的fasta文件,可以使基因组序列也可以是CDS,peptide,protein, 按道理只要是个fasta他都能建立
# -parse_seqids: 对于fasta文件会分析id
# -hash_index: 建立哈希索引
# -dbtype: 数据库类型
  • 本地blast
    我最常用的blastn
blastn -db Dt3rd.fasta -query dtsub.fa -out Dt-result6.xls -evalue 0.00001 -max_target_seqs 1 -num_threads 4 -outfmt "6 qacc sacc evalue length pident "
# -db database
# -query 需要blast的文件
# -out output
# -evalue e值越小越好,最小为0
# -max_target_seqs 最大匹条目
# -num_threads 要使用几个核心
# -outfmt 输出格式

更多关于blast参数和原理参考:
Reference:
blast原理 by 简书 leoinUSA
blast参数 by 博客 牛牛龙

  • for循环对路径下的文件进行批量操作:
    例子: 比如我想做venn,想从差异基因文件中批量提取ID,我一共有17个样本3个组织,如果一个个cp估计自己把自己也整傻了,所以用for循环对路径下的文件批量awk ’{print $1}',然后生成一个个小文件保存起来,然后用TBtools中的venn或者upsetplot拖进去做,速度杠杠的😁
for i in `ls *.xls`; do awk '{print $1}' $i > tmp/$i; done
# 对于路径下.xls的文件,依次将第一列打印出来再保存到 tmp文件下名字还是原来的文件名

当然这只是一个简单的例子,习惯了for后面可能在学习while等。
Oct 23,2019补充

  • while循环批量操作:
    来源于生信技能树Jimmy的教程,这里只粘贴代码,添加注释,具体的例子移步下面reference的链接。

这个脚本有需要注意以下几点:

  1. 需要做一个配置文件,例如这里Jimmy的config包含了3列信息,第一列是样品名字,第二列是第一个文件,第二列是第二个文件,当然,如果你的脚本需要多个文件,你也可以多加几列。
sample file1 file2
SampleA A.file1.txt A.file2.txt
SampleB B.file1.txt B.file2.txt
  1. 这样做有一个好处,就是可以事先在config中管理好你要批处理的文件,提前写好相关信息,在后来output的时候直接用第一列中的变量做标记给结果文件命名。此外,这种方式有助于最后出错的时候纠错,查看配置文件就可以排查出错误。
  2. 在生成config时,可以直接ls 绝对路径/file > filename.txt来保存文件的绝对路径带文件名,这样在操作的时候不会出现路径错误导致报错。
while read id
do
    sample=$(echo $id |cut -d" " -f 1 ) #可以理解成上面for循环中的 i
    file1=$(echo $id |cut -d" " -f 2 ) #这里相当于读入了一个整config文本,那么id对应的就是这个cat config,这里是吧config的第二列提取了出来作为sample这个变量。
    file2=$(echo $id |cut -d" " -f 3 ) 
    echo $sample #echo后每完成一个循环会在shell中打印出该样本的名称。这样就可以看到到哪个样本了,当然,这里在最后加了&,挂在后台,也可以开头打上nohup,最后在nohup里面查看进度。
    echo $file1 
    echo $file2 
~/biosoft/HISAT/current/hisat2 -p 4 --dta -x ~/reference/index/hisat/hg19/genome -1 $file1 -2 $file2 -S $sample.hisat2.hg19.sam 2>$sample.hisat2.hg19.log & # 执行命令
done <$1#和之前提到的一样,$1就是默认的脚本后跟的第一个变量,理解为那个config

另外一个简单的
这里的config只有一列,就是input filename,这里sample相当于把文件名中第一个.前面的保留,比如你的文件都是以sampleName.tissue.xls命名的,那么sample这个变量就是每个文件的sampleName,这里需要根据实际情况自己写。

while read id
do
    file=$(basename $id ) 
    sample=${file%%.*} 
    echo $id $sample 
    nohup samtools sort -@ 4 -o ${sample}.sorted.bam $id & 
done <$1

Reference
Source: 简书,by 生信技能树

  • 没有图形界面(GUI)的Linux下修改配色。

之前自己鼓捣Linux在terminal中自己可以设置颜色,但是用了服务器后发现黑底蓝字看的眼睛快瞎了,再加上我搞了Item2,又加了个骚不拉几的背景图片,更看不清了,一顿百度都是通过 .bashrc下的PS=来修改路径各级颜色的,这不是我想要的,最终终于找到一个教程,简单粗暴好理解:

# 第一步把修改文件夹颜色的文件拷到自己用户的家目录下
cp /etc/DIR_COLORS ~/.dir_colors
# 第二步,修改文件夹颜色的命令
vi ~/.dir_colors ##找到下面的文字
DIR 01;34     # directory
# 将其修改成
DIR 01;33     # directory    
# 即,将深蓝色改成黄色
# 运行bash命令,或者重新登录,即logout然后login

Source: CSDN, by yas_xi

最后编辑于
©著作权归作者所有,转载或内容合作请联系作者
  • 序言:七十年代末,一起剥皮案震惊了整个滨河市,随后出现的几起案子,更是在滨河造成了极大的恐慌,老刑警刘岩,带你破解...
    沈念sama阅读 203,098评论 5 476
  • 序言:滨河连续发生了三起死亡事件,死亡现场离奇诡异,居然都是意外死亡,警方通过查阅死者的电脑和手机,发现死者居然都...
    沈念sama阅读 85,213评论 2 380
  • 文/潘晓璐 我一进店门,熙熙楼的掌柜王于贵愁眉苦脸地迎上来,“玉大人,你说我怎么就摊上这事。” “怎么了?”我有些...
    开封第一讲书人阅读 149,960评论 0 336
  • 文/不坏的土叔 我叫张陵,是天一观的道长。 经常有香客问我,道长,这世上最难降的妖魔是什么? 我笑而不...
    开封第一讲书人阅读 54,519评论 1 273
  • 正文 为了忘掉前任,我火速办了婚礼,结果婚礼上,老公的妹妹穿的比我还像新娘。我一直安慰自己,他们只是感情好,可当我...
    茶点故事阅读 63,512评论 5 364
  • 文/花漫 我一把揭开白布。 她就那样静静地躺着,像睡着了一般。 火红的嫁衣衬着肌肤如雪。 梳的纹丝不乱的头发上,一...
    开封第一讲书人阅读 48,533评论 1 281
  • 那天,我揣着相机与录音,去河边找鬼。 笑死,一个胖子当着我的面吹牛,可吹牛的内容都是我干的。 我是一名探鬼主播,决...
    沈念sama阅读 37,914评论 3 395
  • 文/苍兰香墨 我猛地睁开眼,长吁一口气:“原来是场噩梦啊……” “哼!你这毒妇竟也来了?” 一声冷哼从身侧响起,我...
    开封第一讲书人阅读 36,574评论 0 256
  • 序言:老挝万荣一对情侣失踪,失踪者是张志新(化名)和其女友刘颖,没想到半个月后,有当地人在树林里发现了一具尸体,经...
    沈念sama阅读 40,804评论 1 296
  • 正文 独居荒郊野岭守林人离奇死亡,尸身上长有42处带血的脓包…… 初始之章·张勋 以下内容为张勋视角 年9月15日...
    茶点故事阅读 35,563评论 2 319
  • 正文 我和宋清朗相恋三年,在试婚纱的时候发现自己被绿了。 大学时的朋友给我发了我未婚夫和他白月光在一起吃饭的照片。...
    茶点故事阅读 37,644评论 1 329
  • 序言:一个原本活蹦乱跳的男人离奇死亡,死状恐怖,灵堂内的尸体忽然破棺而出,到底是诈尸还是另有隐情,我是刑警宁泽,带...
    沈念sama阅读 33,350评论 4 318
  • 正文 年R本政府宣布,位于F岛的核电站,受9级特大地震影响,放射性物质发生泄漏。R本人自食恶果不足惜,却给世界环境...
    茶点故事阅读 38,933评论 3 307
  • 文/蒙蒙 一、第九天 我趴在偏房一处隐蔽的房顶上张望。 院中可真热闹,春花似锦、人声如沸。这庄子的主人今日做“春日...
    开封第一讲书人阅读 29,908评论 0 19
  • 文/苍兰香墨 我抬头看了看天上的太阳。三九已至,却和暖如春,着一层夹袄步出监牢的瞬间,已是汗流浃背。 一阵脚步声响...
    开封第一讲书人阅读 31,146评论 1 259
  • 我被黑心中介骗来泰国打工, 没想到刚下飞机就差点儿被人妖公主榨干…… 1. 我叫王不留,地道东北人。 一个月前我还...
    沈念sama阅读 42,847评论 2 349
  • 正文 我出身青楼,却偏偏与公主长得像,于是被迫代替她去往敌国和亲。 传闻我的和亲对象是个残疾皇子,可洞房花烛夜当晚...
    茶点故事阅读 42,361评论 2 342

推荐阅读更多精彩内容