R语言 |一些常用的数据整理的技巧（二）

1. 替换

data$V5[data$V5=='Retroposon']<-"SVA"  #把某列中某值替换成指定值
1

2. 分组计数

#方法一：R中自带的aggregate()函数
group_mean <- aggregate(weight ~ feed, data = df, mean) #这个似乎实现的需求比较简单，只能按照1列进行求和或其他


#方法二：dplyr包中的count()函数
library(dplyr)
selDat<-filter(.data=data,V5=="DNA")
selDat %>% count(V1,V4,V5,name = 'cnt', sort = TRUE) #count()函数只能计数，如果要分组进行其他处理还需要用其他的函数

1
2
3
4
5
6
7
8
9

3. 批量读取目录下的文件

setwd("/home/xxzhang/workplace/project/multiomics/input/V1C-L4IT-AUPs-MIR_HM/")
#install.packages("ggrepel")
library(vroom) #主要使用的package是vroom
library(stringr)
library(ggrepel)
library(reshape2)
filelist <- list.files("./") #需要制定一个目录
n <-length(filelist) 

resultDat<-as.data.frame(matrix(NA,nrow =0,ncol = 2))

for (i in 1:n) #常跟for循环联用，批量快速读取目录下的文件
{
  test<- vroom(filelist[i],id="./",col_names= T,delim="\t")
  t1<-as.character(test[13,2])
  files<-as.character(test[1,1])
  P_value<-strsplit(t1,"\t")[[1]][2]
  ratio<-strsplit(t1,"\t")[[1]][4]
  tissue<-strsplit(files,"_Enhancers.fisher.txt")[[1]][1] #也常用stringr来提取文件名做行名
  values<-cbind(P_value,ratio)
  rownames(values)<-tissue
  resultDat<-rbind(resultDat,values)
}
resultDat$P_value=(-log10(as.numeric(resultDat$P_value)))
resultDat[resultDat$P_value=="Inf","P_value"]=max(resultDat$P_value[resultDat$P_value!=max(resultDat$P_value)])+10

resultDat$label<-as.character(as.data.frame(str_split(rownames(resultDat),"\\."))[2,])
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27

4. cmd指令嵌入perl代码中，实现批量shell

system_call("giggle search -i /home/xxzhang/workplace/project/multiomics/ref/HistomeMarker_index/ -q \"/home/xxzhang/workplace/project/multiomics/Enrichment_cCRE/Enrichment_input-TE_gz/TE-".$ID[0]."-".$sub.".bed.gz\" -s >./giggle_TE_HM/".$ID[0]."-".$sub.".HM.giggle.result");
#在字符串中，用.来连接变量；如果遇到双引号等，为了和字符串本身的双引号区分，要加斜杠\"
sub system_call
{
  my $command=$_[0];
  print "\n\n".$command."\n\n";
  system($command);
}
1
2
3
4
5
6
7
8

5. perl中，使用Getopt::Long模块来获取用户的输入参数

use Getopt::Long;
GetOptions(
            "c=s" =>\$celltype_files,
            "h|help" =>\$help,
);
open (CELL, "< ".$celltype_files) or die "can not open it!"; #$celltype_files即为从用户那边获取的-s后的参数变量
1
2
3
4
5
6

6. PBS系统的示例提交脚本

参考：https://blog.csdn.net/kunxitoothache/article/details/109897918

#PBS -N EP-pair
#PBS -q batch #节点名称，如果占用的内存比较大，需要到gpu或fat节点
#PBS -l walltime=720:00:00 #程序最大运行时间，适用于当自己的任务超过系统常规的任务时，需要自己人为设定；
#PBS -l nodes=1:ppn=16 #规定使用的节点数 也可以指定节点#PBS -l nodes=gpu03:ppn=4
#PBS -o /home/xxzhang/data/Sequencing/scRNA-seq/Zhu_4Region_lateAdulthood/split-fastq/fastq/pbs/tmp/tmp.dlPFC_donor53.out #输出文件地址
#PBS -e /home/xxzhang/data/Sequencing/scRNA-seq/Zhu_4Region_lateAdulthood/split-fastq/fastq/pbs/tmp/tmp.dlPFC_donor53.err #错误文件地址
#PBS -M 2456392738@qq.com #邮件发送的邮箱地址


cd /home/xxzhang/workplace/project/multiomics/E-P-Pair/tmp/result/Part5/
python alignEP.py #注意记得加空格，要不然会出现提交到pbs之后，被终止，什么也没有输出，找不到出错原因#这个坑之前也踩过
1
2
3
4
5
6
7
8
9
10
11

假设上述文件的名称为test.pbs。

qsub test.pbs #把任务命令提交到服务器上
qalter -m ae test.pbs #修改test.pbs任务的邮件发送的属性
qdel test.pbs #删除任务 
1
2
3

7. 使用bioPython的一些模块进行操作

from Bio import pairwise2
from Bio.Seq import Seq 
from Bio.pairwise2 import format_alignment #双序列比对
seq1 = Seq("AGCCTGAGCAACATAGTGAGACCTTGTCTCTACAAAAAATTCAAAACTTAGCCAAGTGTAGTGGTAGGTGCCTTTAGGTGCAGCCACTCAGGAGGCTGAGGTAGGAGGGCTGCTTGAGCCCAGGAGTATGAGGCTGTAGTGAGCTATGATGGTGCCATTGCACTCCAACCTGGACAACAGAGCAAGATCCTGTCTTAAAAAAAGAAAAGAAAA")
seq2 = Seq("GGCCGGGCACAGTGGCTCATGCCTGTAATCCCAGCACTTTGGGAGGCTGAGACGGCAGATCACTGGAGGTCAGGAGTTTGAGACCAGCCTGGCCAACATGGTGAAACCCCGTCTCTACTAAAAATACAAAAATTAGCTGGATGTGGTGATGCTTGCCTGTAATCCTAGCTACTTGGGAGGCTGAGGCAGGGGAATGACTTGAACCTGGGAGGTGGAGATTTTACTCAGCCGAGATGGTGCCATTACACTCCAGCCTGGGTGACAGAGTGAGACTCTGTCTCAAAAA")
alignments = pairwise2.align.localxx(seq1, seq2) #局部比对
for alignment in alignments: 
	print(format_alignment(*alignment))
1
2
3
4
5
6
7
8

8. 提取list的每个第1位

data$TF<-sapply(TF,"[[",1) #提取list TF的每个第1位
1

9. 将数据指定的索引（chr1-chrX）顺序排列----用bedtools

bedtools sort -i A.bed -faidx ../names.txt >B.sorted.bed #names.txt即为指定的顺序 
#之所以这个专门写一节，是因为之前在对bed文件进行排序的时候，走过很多的弯路。用shell自己的sort函数，需要费很大的劲儿才能达到自己想要的结果；
1
2

10. bedtools的slop和flank函数的区别

10.1 slop函数–左右延伸Xbp（包含）

#slop函数和flank函数，在我之前想取gene的promoter和downstream的位置的时候，经常会被混用
#slop函数，是包含已有的bed区域，左右延伸指定的位置
#flank函数，是不包含已有的bed区域，左右延伸指定的位置

cat A.bed #bed文件
chr1 100 200 a1 1 +
chr1 100 200 a2 2 -

cat my.genome #基因组的染色体大小
chr1 1000

bedtools slop  -i A.bed -g my.genome -l 50 -r 80 -s #此处要特别注意，如果按照链进行取的话，链应该在第6列
chr1 50  280 a1 1 +
chr1 20  250 a2 2 -
#这里曾经踩过很严重的坑，因为忽略到strand应该在第6列，所以对负链的处理无效，全部被bedtools当做了正链来处理，因此影响到了下游的一系列的操作；
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15

10.2 flank函数–左右延伸Xbp（不包含）

cat A.bed #bed文件
chr1 100 200
chr1 500 600

cat my.genome #基因组的染色体大小
chr1 1000

bedtools flank -i A.bed -g my.genome -b 5 #both两段各延伸5bp
chr1  95      100
chr1  200     205
chr1  495     500
chr1  600     605

bedtools flank -i A.bed -g my.genome -l 2 -r 3 #left 2bp,right 3bp
chr1  98      100
chr1  200     203
chr1  498     500
chr1  600     603
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18

11. 批量根据某一列的信息拆分大的文件

#!perl
use Getopt::Long;
GetOptions(
            "i=s" =>\$file1,
            "f=s" =>\$file2,
            "e=s" =>\$index,
            "o=s" =>\$type,
            "p=s" =>\$col,
            "h|help" =>\$help,
);
 
print $file1;
print $type;
print $col;

if($help)
{
print
"
usage:
-i         : Hs_repeat_superfamily.txt
-f         : repeat_interest_uniq.sort.bed
-e         : ../split_sort_b
-o         : superfamily
-p         : 5
-h         : usage of this perl file
"
}
##获取用户的输入文件和参数

open (MARK, "< ".$file1) or die "can not open it!"; #读取文档，每一行是我们要根据其提取的特征
while ($line = ){
		print($line);
		chomp($line);
		print($line);
		system_call("mkdir  ./".$type."/".$line); #新建文件夹
		system_call("awk \'\$".$col."==\"".$line."\"\' ".$file2." |sed 's\/\\s\\+\/\\t\/g' >./".$type."/".$line."/".$line.".bed"); #提取文件中，某列为XX的行
}
close(MARK);
sub system_call
{
  my $command=$_[0];
  print "\n\n".$command."\n\n";
  system($command);
}

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46

12. 使用Homer做motif富集

(1)构建homer的输入文件

awk -v FS="," -v OFS="\t" '{if($13~"TRUE") print "peak"NR,$5,$6,$7,$8 }' AU-L1_in_dlPFC_L2-3_IT.csv >Hs_AU-L1.bed
awk -v FS="," -v OFS="\t" '{if($13~"FALSE"||$13~"TRUE") print "peak"NR,$5,$6,$7,$8 }' AU-L1_in_dlPFC_L2-3_IT.csv >AU-L1.bed
1
2

主要要注意的是“peak”NR的应用；

(2) 使用homer进行motif富集

findMotifsGenome.pl /home/xxzhang/workplace/project/multiomics/input/AUPs/Homer_input/dlPFC-L2-3-IT-AUPs-Homer.bed  hg38  /home/xxzhang/workplace/project/multiomics/output/HomerTF/  -size given
1

13. 去除全是零的行

X[which(rowSums(X) > 0),]
1

14. PCA 分析

filter_data2<-filter_data[!apply(filter_data,1,var)%in%c(0,NA),]
pca <- prcomp(t(filter_data2), center=TRUE, scale=TRUE) #注意这里需要转置

pcaPlotDat <- cbind(filter_phe, pca$x[,1:2])
ggplot(pcaPlotDat) + geom_point(aes(x=PC1, y=PC2, shape=Species, color=factor(Predicted.period)), size=3)
1
2
3
4
5