*代码笔记

作者: 曲凉不见 | 来源:发表于2019-12-31 16:52 被阅读0次

一、服务器用户端

NT/NR库：
ftp://ftp.ncbi.nlm.nih.gov/blast/db/FASTA/

makeblastdb -in nt -parse_seqids -hash_index -dbtype nucl -logfile nt_logfile
makeblastdb -in nr -parse_seqids -hash_index -dbtype prot -logfile nr_logfile

二、Unix/Linux命令

查看硬盘分区及使用情况
df -lh

1. vi下:

当前行行首：0（数字0）
当前行行尾：$
文件尾：G
当前屏首行：H
当前屏末行：L
显示行号：":set number"

2. 四则运算

echo "xxx" | bc -l #双引号里面写运算

3. echo命令

echo打出tab
echo -e "A\tB\tC" #echo完了自带换行符

wget下载限制
wget -c --limit-rate=3M [http://hdahdoad/file](http://hdahdoad/file) #断了能接上，且最大速度不超过3M

4. awk命令

输出多列：
cat file | awk -F "\t" '{print $1,$2,$3}' OFS="\t"

删除重复行：
awk '{if(!seen[$0]++){print $0;}}' file

删除最后一列：
awk '{print $NF}' file #$NR是最后一行

删除第一行：
awk '{$1=" ";print $0}' file

求和/求均值：
awk '{sum+=$1}END{print "sum=",sum}' file
awk '{if($1>=5)sum+=$1}END{print sum}' file
awk '{ sum += $7; } END { print "sum = " sum; print "average = " sum/NR }'

这是什么：
awk '{a[NR]=length($0)}END{for(i in a) if(i%4==2 && a[i]!=a[i+2])print i"\t"a[i]"\t"a[i+2]}'

CSV转TXT
awk 'BEGIN{FS=",";OFS="\t"}{$1=$1}1' file

FASTQ文件,每个reads输出序列行
awk '{if(NR%4==2)print $0}' test.fastq #NR是行，length($0)是长度

计算第二列之和
awk '{sum+=$2;count+=1}END{print "SUM:"sum"\nAVG:"sum/count"\nCOUNT:"count}' file

去掉fasta中多余的换行符
awk '!/^>/{printf "%s",$0;n="\n"}/^>/{print n $0;n=" "}END{printf "%s",n}' file #小写s

awk与正则表达式的结合
取第三列以chr21开头的reads打印出来，用\~表示匹配，\~!等同\!~
awk '{if($3~/^chr21/)print $0}' sample.all.1.sam | less -S

/REG/为正则表达式，可以将$0中满足条件的记录送到action处理。
awk '/REG/{action}'

以Permission denied结尾的行输出
awk '{if($0~/Permission denied$/)print $0}' file

awk中使用if、elsif和else
awk -F '\t' '{if($5=="-")print $1,$2,$3,$4,$6,$6,$7,$8,$9;else print $1,$2,$3,$4,$5,$6,$7,$8,$9}' OFS="\t" NC_000962.ptt > NC_000962_re.ptt

统计出现次数
awk '{a[$3]++} END{for(i in a){print i,a[i] | "sort -r -n -k 2"}}' file
这是统计文件第三列的信息，重复出现的次数，输出第三列，以及次数，次数在结果的第二列，因此对结果的第二列进行排序输出

5. tr命令（替换）

去掉大写C：
tr -d "C" file

A换成B：
tr -s "A" "B" file

6. sed命令

替换（去掉）分号：
sed -i 's/;//g' file #-i是不会生成新的文件，而是生成是文件直接代替file

取出文件中的第n行
sed -n '5,10p' file > output  #取出5-10行,不改变源文件

去掉含有某些字符的行
sed -i '/string/d' file

删除第n行
sed 'nd' file

从第三行开始，每隔一行删一行
sed 3~2d file

删除第4-8行
sed 4,8d file > output (前面加-i 就直接修改文件本身了)

删除最后一行
sed '$'d file

匹配行删除
sed /string/d file

删除从匹配行到末尾
sed '/string/,$d' file

删除从匹配行到之后两行
sed '/string/,+2d' file

删除空行
sed '/^$/d' file

在文件中增加一行
sed -i 'echo\"2\";/a\echo \"5\";' file #在echo "2";下面增加一行是echo "5";

替换大小写
sed 's/[A-Z]/\l&/g' file > file_cout #变小写
sed 's/[a-z]/\L&/g' file > file_cout #变大写

在第一行前面插入一行
sed -i '1iaaa\tbbb\tccc\tddd' file #在file的第一行增加aaa\tbbb\tccc\tddd

7. sort命令

按第x列排列（从小到大）：
sort -k -x -n file

8. grep命令

精确抓取：
grep "\<AAA\>" file #只写一半就是半边精确匹配的抓取

范围抓取：
grep "AAA[12]" file #抓AAA1或AAA2

grep同时抓取多个关键词
grep -E "AA|BB|CC" file

9. find命令

find /path/../ -type f -name "*.pl"

10. touch命令

touch创建文件或修改文件时间

选项
-a 只更新访问时间，不改变修改时间
-c 不创建不存在的文件
-m 只更新修改时间，不改变访问时间
-r file 使用文件file的时间更新文件的时间
-t 将时间修改为参数指定的日期,如：07081556代表7月8号15点56分

Example:
touch -t 201601291200.00 reads.split.050.2.fastq

11. du/df命令

查看当前目录下所有文件和文件夹大小（不单独列出单个文件大小）
du -h

查看目标文件夹大小
du -sh /path/

12. 从github下载软件安装包

假设一个软件的github网页为：https://github.com/xavierdidelot/ClonalFrameML

git clone [https://github.com/xavierdidelot/ClonalFrameML](https://github.com/xavierdidelot/ClonalFrameML)

git clone [https://github.com/xavierdidelot/ClonalFrameML.git](https://github.com/xavierdidelot/ClonalFrameML.git)

13. nohup

输出log文件到指定文件中
nohup sh xxx.sh >> log &

14. md5sum

输出MD5
md5sum file

检测MD5
md5sum -c md5.check.txt

md5.check.txt:
12345 file1
56789 file2
中间是两个空格隔开

15. wget

断点续传
wget -c ftp://xxx.fq.gz

下载文件夹
wget -c -r -np -nH ftp://xxx
-r : 遍历所有子目录
-np : 不到上一层子目录去
-nH : 不要将文件保存到主机名文件夹
-R index.html : 不下载 index.html 文件

三、Shell命令

高亮grep结果，在~/.bashrc里加一行

Alias grep='GREP_COLOR="1;33;40" LANG=C grep--color=auto'
source ~/.bashrc

循环
file.txt里面是文件名，把这些文件的第二列都拿出来放在zzz里

for filename in $(cat file.txt)
do
  samtools view 1-mapping/$filename | cut -f2 >>zzz
done

四、Python

1. 安装python模块

安装到root下

python -m pip install SomePackage

安装到用户自己目录下

python -m pip install SomePackage --user

检查某模块的版本

>>> import pkg_resources
>>> pkg_resources.get_distribution("SomePackage").version

安装本地模块

pip install ./downloads/SomePackage-1.0.4.tar.gz

或卸载

pip uninstall package-name

五、Perl语言

1. Last, next, break

Perl中last, next, break的用法：
Next是跳到下一次循环i+1，last是跳出循环。

2.去掉字符左右的空白

s/(^s+|s+$)//g

3.去掉哈希中的一个键值对

delete $hash{$key}

4.操作数组

4.1 从数组右端操作：

取值：pop
增值：push @array, $a;

4.2 从数组左端操作

取值：shift
增值：unshift

4.3 foreach遍历array中的值

foreach $rock(@rocks){
 $rock="\t$rock";
 $rock="$rock\n";
}

4.4 数组排序：sort不会改变array本身，所以需要保存。

@sorted=sort(@rocks);
@back=reverse sort @rocks;
@rocks=sort @rocks; #重新存回也可以

4.5 数字排序会变成ASCII码来，1开头的会在3开头的前面，解决方法：

@num=sort {$a<=>$b} @array; #从小到大
@num=sort {$b<=>$a} @array; #从大到小

4.6 数组值赋给哈希值

$hash{$key}=[@array];
print "@{$hash{$key}}\n";

4.7 数组去重

my %count;
my @new_array = grep { ++$count{ $_ } < 2; } @array;

5.数组与字符串的转换

用tab键隔开的字符串
$chr=join("\t",@array)

这样就是直接变成array
@array=split /\t/

6. 数组去重

my @array = ( 'a', 'b', 'c', 'a', 'd', 1, 2, 5, 1, 5 );
my %saw;
@saw{ @array } = ( );
my @uniq_array = sort keys %saw;

7. 按照键/值的顺序输出哈希

7.1. 按照key 排序
7.1.1 按照数字排序

foreach my $item (sort {$a<=>$b} keys %hash){
 print "$item == > $hash{$item}","/n";
}

7.1.2 按照ASCII排序

foreach my $item (sort {$a cmp $b} keys %hash){
 print "$item == > $hash{$item}","/n";
}

7.2.按照value排序
7.2.1 按照数值排序（大到小，小到大只需要将 $a 和$ b 位置调换即可）

foreach my $key ( sort { $hash{$a} <=> $hash{$b} } keys %hash ) {
 my $value = $hash{$key};
 #do something with ($key, $value)
}

7.2.2 按照ACSII排序

foreach my $key ( sort { $hash{$a} cmp $hash{$b} } keys %hash ) {
 my $value = $hash{$key};
 # do something with ($key, $value)
}

8. Substr获取字符串

offset代表起始字符的位置,length代表引用的字符串长度，如果省略length则代表从起始值到字符串的最后一个字符长度。而offset如果是负值的话，就会从字符串右边开始指定字符。

substr($string,offset,length)
$s=substr("perl5",2,2);        这时$s="rl"
$s=substr("perl5",2);          这时$s="rl5"
$s=substr("perl5",-2,2);       这时$s="l5"

9. Splice去除数组中的值

splice(@array,skipped,length,@newarray)
skipped指，在去掉某个值之前跳过几个值
length指，去掉几个值

10. 正则表达式计算一段字符中，匹配到短字符的次数

计算字符串 $var 中出现的A字符，用到了正则匹配"s///g"
my $var = 'TCTCATGTGAAAAACTATATCAATAATATAAAAACA';
my $count = ($var =~ s/A/A/g);
print $count;

11. 正则表达式计算一段字符中，获取匹配指定字符的部分，返回其位置

while ($sequence=~m/(N+)/g) {# 匹配一个N或以上的字符
    my $len=length($1);# 返回这段匹配的长度
    my $end=pos($sequence);#用 pos函数返回该匹配的终止位置
    my $start=$end-$len+1;# 计算出起始位置
    print OUT "$key\tN\t$start\t$end\t$len\n";# 输出结果
}

12. 一个字符串，替换字符，写成函数

sub write_string {
    my ($sequence,$ref_seq,$start,$end)=@_;
    my $len_before=$start-1;
    my $string_before=substr($sequence,0,$len_before);

    my $offset=$start-1;
    my $len_middle=$end-$start+1;
    my $string_mid=substr($ref_seq,$offset,$len_middle);

    my $offset_after=$end;
    my $len_after=length($sequence)-$end;
    my $string_after=substr($sequence,$offset_after,$len_after);

    my $new_string=join("",$string_before,$string_mid,$string_after);
    return $new_string;
}
!使用的时候
my $new_string=write_string($sequence,$ref_sequence,$start,$end);

六、R语言

1. R更新

install.packages("installr")
library(installr)
然后在工具栏选择更新R

改变工作目录：

setwd("\path\...")

查看数据类型

mode(x)

2. apply函数家族

http://blog.fens.me/r-apply/
apply函数族是R语言中数据处理的一组核心函数，通过使用apply函数，我们可以实现对数据的循环、分组、过滤、类型控制等操作。

3. R筛选值

只保留第二列值大于0的行
data=read.table("6-RankSum_11492_1_methy_5k.txt")
data=data[data[,2]>0,]

4. apply函数

apply函数本身就是解决数据循环处理的问题，为了面向不同的数据类型，不同的返回值，apply函数组成了一个函数族，包括了8个功能类似的函数。这其中有些函数很相似，有些也不是太一样的。

apply函数是最常用的代替for循环的函数。apply函数可以对矩阵、数据框、数组(二维、多维)，按行或列进行循环计算，对子元素进行迭代，并把子元素以参数传递的形式给自定义的FUN函数中，并以返回计算结果。

apply(X, MARGIN, FUN, ...)
X:数组、矩阵、数据框
MARGIN: 按行计算或按按列计算，1表示按行，2表示按列
FUN: 自定义的调用函数
…: 更多参数，可选

比如，对一个矩阵的每一行求和，下面就要用到apply做循环了。
> x<-matrix(1:12,ncol=3)
> apply(x,1,sum)
[1] 15 18 21 24

下面计算一个稍微复杂点的例子，按行循环，让数据框的x1列加1，并计算出x1,x2列的均值。
# 生成data.frame
> x <- cbind(x1 = 3, x2 = c(4:1, 2:5)); x
 x1 x2
[1,] 3 4
[2,] 3 3
[3,] 3 2
[4,] 3 1
[5,] 3 2
[6,] 3 3
[7,] 3 4
[8,] 3 5

> data.frame(x1=x[,1]+1,x2=rowMeans(x))
 x1 x2
1 4 3.5
2 4 3.0
3 4 2.5
4 4 2.0
5 4 2.5
6 4 3.0
7 4 3.5
8 4 4.0

4. ggplot

col=sample是连续值，col=factor(sample)是不连续值

4.1 RColorBrewer(??RColorBrewer)

library(RColorBrewer)
mypalette<-rev(brewer.pal(9,"Greens"))#allowed maximum for palette Greens is 9
pheatmap(data.rev,col=mypalette,cluster_col=F,cluster_row=F)

Creates nice looking color palettes especially for thematic maps

brewer.pal(n, name)
display.brewer.pal(n, name)

n	Number of different colors in the palette, minimum 3, maximum depending on palette
name	A palette name from the lists below
type	One of the string "div", "qual", "seq", or "all"
select	A list of names of existing palettes
exact.n	If TRUE, only display palettes with a color number given by n
colorblindFriendly	if TRUE, display only colorblind friendly palettes

一般用法可以是：

library(RColorBrewer)
myPalette <- colorRampPalette(rev(brewer.pal(11, "RdBu")))
p=ggplot(sample,aes(x = X2, y = X1, fill = value))
p=p + geom_tile()
p=p + scale_fill_gradientn(colours = myPalette(100))

4.2 Facet(分面)
虽然我们前面说过ggplot2分面最终的效果是一页多图，但跟通常所说的在 “一个页面中绘制多个图形”还是有区别的。ggplot2分面体现的是数据分组，同一页面中的多个小图是完全相同的类型。真正意义的“一页多图”在 ggplot2中需要通过其他方法实现。ggplot2的分面有两种方式，分别使用facet_wrap或facet_grid 函数。

facet_grid(.~sample)     横着排
facet_grid(sample~.)     竖着排

缠绕分面 facet_wrap
facet_warp 即“缠绕分面”，对数据分类只能应用一个标准，不同组数据获得的小形按从左到右从上到下的“缠绕”顺序进行排列：

facet_wrap(facets, nrow = NULL, ncol = NULL, scales = "fixed", shrink = TRUE, as.table = TRUE, drop = TRUE)

facets	分面参数如 ~cut，表示用 cut 变量进行数据分类
nrow	绘制图形的行数
ncol	绘制图形的列数，一般nrow/ncol只设定一个即可
scales	坐标刻度的范围，可以设定四种类型。fixed 表示所有小图均使用统一坐标范围；free表示每个小图按照各自数据范围自由调整坐标刻度范围；free_x为自由调整x轴刻度范围；free_y为自由调整y轴刻度范围。
shrinks	也和坐标轴刻度有关，如果为TRUE（默认值）则按统计后的数据调整刻度范围，否则按统计前的数据设定坐标。
as.table	和小图排列顺序有关的选项。如果为TRUE（默认）则按表格方式排列，即最大值（指分组level值）排在表格最后即右下角，否则排在左上角。
drop	是否丢弃没有数据的分组，如果为TRUE（默认），则空数据组不绘图。

格网分面 facet_grid
格网分面可以应用多个标准对数据进行分组。
color~cut 对数据的分组和小图排列有决定作用，波浪号前为小图分行标准，后面为分列标准。facet_grid 的完整用法为：

facet_grid(facets, margins = FALSE, scales = "fixed", space = "fixed", shrink = TRUE, labeller = "label_value", as.table = TRUE, drop = TRUE)

和facet_wrap比较，除不用设置ncol和nrow外（facets公式已经包含）外还有几个参数不同：

Margins	这不是设定图形边界的参数。它是指用于分面的包含每个变量元素所有数据的数据组
space	这个参数要配合scales使用，如果为fixed（默认），所有小图的大小都一样，如果为

free/free_x/free_y，小图的大小将按照坐标轴的跨度比例进行设置。
labeller|这是设定小图标签的，facet_grid的函数说明档讲得比较明白，参考之。或许会在后面介绍。

4.3 一个画布中画多个ggplot图

library(easyGgplot2)
ggplot2.multiplot(p1,p2,p3,p4, cols=2)

5. R中获取一个data的大小

dim(data) #行列, ncol(data) #列, nrow(data) #行

Melt和cast
melt是使矩阵变少一维
cast是melt的反转，例如

miRNA	Gene	Estimate
hsa-miRNA-1	AA	-0.91
hsa-miRNA-2	BB	-0.99
hsa-miRNA-1	AA	-0.97
hsa-miRNA-2	BB	-0.95

acast(data,Gene~miRNA,mean)变成

Gene	hsa-miRNA-1	hsa-miRNA-2
AA	-0.91	-0.97
BB	-0.99	-0.95

当不仅有miRNA，Gene的时候，Estimate会有不止一个值，mean就是取平均。

6. R中用0代替NA

data.c[is.na(data.c)]=0

7. R中取特定的行

而利用subset()函数进行访问和选取数据框的数据更为灵活，subset函数将满足条件的向量、矩阵和数据框按子集的方式返回。

Subset函数的三种应用方式：

subset(x, subset, ...)                     
subset(x, subset, select, drop = FALSE, ...)    ##对于矩阵
subset(x, subset, select, drop = FALSE, ...)    ##对于数据框
x是对象，subset是保留元素或者行列的逻辑表达式，对于缺失值用NA代替。
Select 是选取的范围，应小于x。 x<-data.frame(matrix(1:30,nrow=5,byrow=T))
rownames(x)=c("one","two","three","four","five")
colnames(x)=c("a","b","c","d","e","f")
x
new<-subset(x,a>=14,select=a:f)
new                            ## 从a到f列选取a>14的行