注册    登录
专栏名称: ShawnMagic
什么时候可以毕业ಥ_ಥ
今天看啥
微信公众号rss订阅, 微信rss, 稳定的RSS源
目录
相关文章推荐
今天看啥  ›  专栏  ›  ShawnMagic

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

ShawnMagic  · 简书  ·  · 2019-04-12 19:22

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:
😹😹，暂时忘了，以后补上

  • 批量删除环境变量中的数据： 做完数据清洗往往一些很大的原始数据还占据在环境变量中，所以想要删除，清空的命令 rm(list = ls()) ，但是需要删除一些特定的含有共同字符的怎么办？刚开始想到通配符* rm(lncRNA.*) ，当然这个肯定报错，然后在高人指点下完成了，在此感谢BMK原老师🌹 
rm(list=ls(pattern="lncRNA."))# 比如我之前导入的rawdata是lncRNA的FPKM，批量删除那么就让pattern="lncRNA"

An other solution rm(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 ，造成在与其他 database ID不统一，存在 .1 , .2 等 isoform 的标记，其实这里需要做的就是删除结尾的两个字符就可以，这里用 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




原文地址：访问原文地址
快照地址： 访问文章快照
推荐文章
第一财经资讯  ·  东航空难搜救最新进展：未放弃找寻幸存者，机组没有挂出任何紧急代码
2 年前
绿茶书情  ·  2020年度好书榜速写 | 茶画会
3 年前
环境经济观察  ·  《绿色金融理论与实务》最新课程（17-18）
4 年前
ioncology  ·  BC青年团︱李江涛：《浸润性乳腺癌患者引发的思考》
5 年前
吉他谱  ·  蔡老师吉他课堂第四十三讲
6 年前
关于   移动版
今天看啥 - 微信公众号rss订阅, 微信rss, 稳定的RSS源
© 2024 ~ 沪ICP备11025650号