ViewBS:DNA甲基化数据的可视化

官网提供了三种安装方式:conda、docker和正常一步步安装,前两种较为简单,能解决较多的依赖关系,下边介绍一步步安装(https://github.com/xie186/ViewBS)

#安装htslib
git clone https://github.com/samtools/htslib.git
git submodule update --init --recursive
autoreconf
./configure
make
sudo make insatll
#安装ViewBS
wget -c https://github.com/xie186/ViewBS/archive/refs/tags/v0.1.10.tar.gz
tar -zxvf v0.1.10.tar.gz
cd ViewBS-0.1.10
/opt/biosoft/ViewBS-0.1.10/ext_tools/cpanm --local-lib=~/perl5 local::lib && eval $(perl -I ~/perl5/lib/perl5/ -Mlocal::lib)
/opt/biosoft/ViewBS-0.1.10/ext_tools/cpanm Getopt::Long::Subcommand
/opt/biosoft/ViewBS-0.1.10/ext_tools/cpanm Getopt::Long (> 2.50)
/opt/biosoft/ViewBS-0.1.10/ext_tools/cpanm --force Bio::DB::HTS::Tabix

ViewBS主要以Bismark的bismark_methylation_extractor输出结果为输入,同时根据不同的目的需要准备全基因组DNA序列的fasta文件,所有基因的bed文件,转座子bed文件,差异基因的bed文件等。

#数据准备
samtools faidx genome.fasta
bgzip ../A.1_bismark_bt2_pe.deduplicated.CX_report.txt ./ #Bismark结果需要bgzip压缩
tabix -p vcf A.1_bismark_bt2_pe.deduplicated.CX_report.txt.gz #生成tbi结尾的index文件
#分析流程
ViewBS MethCoverage --reference /home/wuchangsong/gc_genome/17.Geta/0.initial_data/genome.fasta --sample A.1_bismark_bt2_pe.deduplicated.CX_report.txt.gz,Control-1 --sample B.1_bismark_bt2_pe.deduplicated.CX_report.txt.gz,Control-2 --sample C.1_bismark_bt2_pe.deduplicated.CX_report.txt.gz,Control-3 --sample D.1_bismark_bt2_pe.deduplicated.CX_report.txt.gz,Single-1 --sample E.1_bismark_bt2_pe.deduplicated.CX_report.txt.gz,Single-2 --sample F.1_bismark_bt2_pe.deduplicated.CX_report.txt.gz,Single-3 --sample G.1_bismark_bt2_pe.deduplicated.CX_report.txt.gz,Multiple-1 --sample H.1_bismark_bt2_pe.deduplicated.CX_report.txt.gz,Multiple-2 --sample I.1_bismark_bt2_pe.deduplicated.CX_report.txt.gz,Multiple-3 --outdir MethCoverage --prefix BS_seq_allsam
ViewBS MethGeno --genomeLength /home/wuchangsong/gc_genome/17.Geta/0.initial_data/genome.fasta.fai --sample A.1_bismark_bt2_pe.deduplicated.CX_report.txt.gz,Control-1 --sample B.1_bismark_bt2_pe.deduplicated.CX_report.txt.gz,Control-2 --sample C.1_bismark_bt2_pe.deduplicated.CX_report.txt.gz,Control-3 --sample D.1_bismark_bt2_pe.deduplicated.CX_report.txt.gz,Single-1 --sample E.1_bismark_bt2_pe.deduplicated.CX_report.txt.gz,Single-2 --sample F.1_bismark_bt2_pe.deduplicated.CX_report.txt.gz,Single-3 --sample G.1_bismark_bt2_pe.deduplicated.CX_report.txt.gz,Multiple-1 --sample H.1_bismark_bt2_pe.deduplicated.CX_report.txt.gz,Multiple-2 --sample I.1_bismark_bt2_pe.deduplicated.CX_report.txt.gz,Multiple-3 --prefix BS_geno_sample --context CG --outdir MethGeno_100k --minLength 1000000 --win 100000 --step 100000
ViewBS MethOverRegion --region repeats.bed --sample A.1_bismark_bt2_pe.deduplicated.CX_report.txt.gz,Control-1 --sample B.1_bismark_bt2_pe.deduplicated.CX_report.txt.gz,Control-2 --sample C.1_bismark_bt2_pe.deduplicated.CX_report.txt.gz,Control-3 --sample D.1_bismark_bt2_pe.deduplicated.CX_report.txt.gz,Single-1 --sample E.1_bismark_bt2_pe.deduplicated.CX_report.txt.gz,Single-2 --sample F.1_bismark_bt2_pe.deduplicated.CX_report.txt.gz,Single-3 --sample G.1_bismark_bt2_pe.deduplicated.CX_report.txt.gz,Multiple-1 --sample H.1_bismark_bt2_pe.deduplicated.CX_report.txt.gz,Multiple-2 --sample I.1_bismark_bt2_pe.deduplicated.CX_report.txt.gz,Multiple-3 --prefix bis_gene_all_sample --context CG --outdir MethOverRegion

详细流程参考官网

awk命令将fastq文件转换成fasta

awk '{if(NR%4==1||NR%4==2){print $0}}'  test.fq | sed 's/^@/>/g' > myfile.fasta

 有时候,比如说我们需要做一个clustalw,我们需要将多行的fasta文件(multiple line fasta)格式文件,转换成单行的fasta文件(single line fasta)序列文件,怎么办呢。 巧用awk + printf一行搞定。
awk '/^>/ { print n $0;}  !/^>/ {printf "%s", $0, n="\n"}  END {print ""}'  test.fa

 

sed,awk

sed工具

sed本身也是一个管道命令,可以分析standard input的,而且sed还可以将数据进行替换、删除、新增、选取特定行等的功能。
用法:sed [-nefr] 动作

参数:
-n:使用安静模式。在一般sed的用法中,所有来自stdin的数据一般都会被列出到屏幕上。但如果加上-n参数后,则只有经过sed特殊处理的那一行(或者操作)才会被列出。
-e:直接在命令行模式上进行sed的动作编辑。
-f:直接将sed的动作写在一个文件内,-f filename则可以执行filename内的sed动作。
-r:sed的动作支持的是扩展型正则表达式的语法(默认是基础正则表达式语法)。
-i:直接修改读取的文件内容,而不是屏幕输出。
动作说明:[n1[,n2]]function
n1,n2:不见得会存在,一般代表选择进行动作的行数,举例来说,如果我的动作是需要在10到20行之间进行的,则‘10,20[动作行为]’
function有下面这些参数:
a:新增,a的后面可以接字符串,而这些字符串会在新的一行出现(目前的下一行);
c:替换,c的后面可以接字符串,这些字符串可以替换n1,n2之间的行;
d:删除,因为是删除,所以d后面通常不接任何参数;
i:插入,i的后面可以接字符串,而这些字符串会在新的一行出现(目前的上一行);
p:打印,将某个选择的数据打印出来。通常p会与参数sed -n一起运行;
s:替换,可以直接进行替换工作。通常这个s的动作可以搭配正则表达式。

例一:删除第2到5行

注意sed后面接的动作,务必以”两个单引号括住!

如果体型变换一下,删除第三到最后一行,则是‘nl /etc/passwd | sed ‘3,$d’’,这个$代表最后一行。

例二:将第2~5行的内容替换成为‘No 2-5 number’

sed 另一个强大的用法,部分数据的查找并替换:sed ‘s/要被替换的字符串/新的字符串/g’

awk:好用的数据处理工具

相比于sed常常作用于一整行的处理,awk则比较倾向于将一行分成数个‘字段’来处理。因此。awk适合小型的数据处理
用法:awk ‘条件类型 1{动作 1} 条件类型 2{动作 2} …’ filename
awk的处理流程:

  1. 读入第一行,并将第一行的数据填入$0,$1,$2等变量当中;
  2. 依据条件类型的限制,判断是否需要进行后面的动作;
  3. 做完所有的动作与条件类型;
  4. 若还有后续的‘行’的数据,则重复上面1-3的步骤,直到所有的数据都读完为止。

变量名字及意义:

  • NF:每一行($0)拥有的字段总数
  • NR:目前awk所处理的是‘第几行’数据
  • FS:目前的分割字符,默认是空格

例:awk的逻辑运算

数据:
处理:

  • 第一行只是说明,所以第一行不要进行加总(NR==1时处理)
  • 第二行以后就会有加总的情况(NR>=2以后处理)
  • 第三行将处理后的结果格式输出
  • 所有的awk动作,即在{}内的动作,如果有需要多个命令辅助时,可利用分号“;”间隔,或者直接以enter键来隔开每个命令
  • 逻辑运算中,如果是“等于”的情况,务必使用两个等号“==”
  • 与bash、shell的变量不同,在awk中,变量可以直接使用,不需要加上$符号

 

以上内容摘自《鸟哥的私房菜》第三版 基础学习篇

排序命令:sort,wc,uniq

sort

用法:sort [-fbMnrtuk] [file or stdin]

参数:
-f:忽略大小写
-b:忽略最前面的空格部分
-M:以月份名字来排序,如JAN,DEC等的排序方法
-n:使用‘纯数字’排序(默认是以文字类型来排序的)
-r:反向排序
-u:就是uniq,相同的数据中,仅出现一行代表
-t:分隔符,默认是用tab键来分割
-k:以那个区间(field)来进行排序

wc

用法:wc [-lwm]

参数:
-l:仅列出行
-w:仅列出多少字
-m:多少字符

uniq

用法:uniq [-ic]

参数:
-i:忽略大小写
-c:进行计数

摘自:《鸟哥的私房菜》第三版 基础学习篇

选取命令:cut,grep

cut

cut -d '分隔字符' -f fields
cut -c 字符范围   <==用于排列整齐的信息

grep

grep [-acinv] [--color=auto] '查找字符串' filename
参数:
-a: 将binary文件以test文件的方式查找数据;
-c: 计算找到‘查找字符串’的次数;
-i: 忽略大小写的不同;
-n: 输出行号;
-v: 反向选择,即显示出没有‘查找字符串’ 内容的那一行;
--color=auto: 可以将找到的关键字部分加上颜色显示;
高级用法:grep [-A] [-B] [--color=auto] '搜寻字符串' filename
参数:
-A:后面可接数字,为after的意思,除了列出该行外,后续的n行也列出来;
-B:后面可接数字,为befer的意思,除了列出该行外,前面的n行也列出来

摘自:《鸟哥的私房菜》第三版 基础学习篇

通配符与特殊符号

通配符是bash操作环境中一个非常有用的功能,利用它我们处理数据就更加方便。

*:代表0个到无穷多个任意字符
?: 代表一定有一个任意字符
[]: 同样代表一定有一个在中括号内的字符(非任意字符)
[-]: 若有减号在中括号内时,代表『在编码顺序内的所有字符』。例如 [0-9] 代表 0 到 9 之间的所有数字,因为数字的语系编码是连续的
[^]: 若中括号内的第一个字符为指数符号 (^) ,那表示『反向选择』,例如 [^abc] 代表 一定有一个字符,只要是非 a, b, c 的其他字符就接受的意思

特殊字符

#: 注释,这个最常用在script中,视为说明.其后的数据均不执行
\: 转义符号,将特殊字符或通配符还原成一般字符
|: 分隔两个管线命令的界定
;: 连续性命令的界定(注意,与管线命令并不相同)
~: 用户的主文件夹
$: 使用变量前导符
&: 将指令变成在背景下工作
!: 逻辑运算中的“非”
/: 路径分隔符号
>,>>: 数据流重定向,输出导向,代表替换和累加
<,<<: 数据流重定向,输入导向
'': 单引号,不具有变量置换的功能
"": 具有变量置换的功能
``: 两个“`”中间为可以先执行的指令
(): 中间为子shell的起始与结束
{}: 中间为命令块的组合

以上为bash环境中常见的特殊字符,理论上,我们的文件名尽量不要使用到上述的字符。

摘自:《鸟哥的私房菜》第三版 基础学习篇

shell的变量功能(二)

变量键盘读取,数组与声明:read,array,declare

1. read

读取来自键盘输入的变量,常被用在 shell script 的撰写当中。

用法:

2. declare / typeset

declare 或 typeset 是一样的功能,就是在声明变量的类型。如果使用 declare 后面并没有接任何参数,那么 bash 就会主动的将所有的变量名称与内容通通叫出来, 就好像使用 set 一样。

用法:declare [-aixr] variable

参数:

-a : 将后面名为 variable 的变量定义成为数组 (array) 类型
-i : 将后面名为 variable 的变量定义成为整数数字 (integer) 类型
-x : 用法与 export 一样, 就是将后面的 variable 变成环境变量;
-r : 将变量设置成为 readonly 类型, 该变量不可被更改内容, 也不能重设

3. array

更多详情请参照《鸟哥的私房菜》第三版 第十一章

与文件系统及程序的限制关系:ulimit

我们可以通过bash 限制用户的某些系统资源,包括打开的文件数量、可以使用的cpu 时间、可以使用的内存总量等。如何设置?用ulimit吧

ulimit [-SHacdfltu] [配额]

-H :hard limit,严格的设置,必定不能超过这个设置的数值

-S :soft limit,警告的设置,可以超过这个设定值,但是若超过则有警告信息。 在设置上,通常soft 会比hard 小,举例来说,soft 可设置为80而hard 设置为100,那么你可以使用到90,但介于80~100之间时,系统会有警告信息通知你。

-a :后面不接任何参数,可列出所有的限制额度;

-c :当某些进程发生错误时,系统可能会将该进程在内存中的信息写成文件(排错用),这种文件就被称为内核文件,此为限制每个文件内核文件的最大容量。

-f :此shell 可以创建的最大文件容量(一般可能设置为2GB)单位为kb

-d :进程可使用的最大断裂内存(segment )容量

-l :可用于锁定(lock)的内存容量

-t :可使用的最大cpu 时间(单位为秒)

-u :单一用户可以使用的最大进程(process)数量;

Ulimit –a

Ulimit –f 10240

Ulimit –a

一般用户使用这个命令的时候只能减少设置的值,不能增加设置的值

 

版权声明:本文为博主原创文章,未经博主允许不得转载。

 

shell的变量功能(一)

什么是变量?

变量是一个存储位置和一个关联的符号名字,这个存储位置包含了一些已知或未知的量或者信息,即值。在C语言里,变量是如下三方面的统一体:

  1. 名字(运行时会变成数字化的名字,内存地址)
  2. 存储位置(某一位置开始的一定大小的存储空间)
  3. 该存储位置里内容的解释方式(即类型,整数、浮点数还是字符串?)

任意一部分单独都不是变量。当我们给一个变量a赋值另一个值时,改变的是a对应的存储位置里的内容,赋值前后是同一个a,因为1、2、3都没有变。

变量的显示

变量的设置规则

  1. 变量与变量内容以一个等号“=”来连接;
  2. 等号两边不能直接接空格符;
  3. shell中所有变量都定义为字符串,且变量名称只能是英文字母与数字,但是数字不能是开头字符;
  4. 若有空格符可以使用双引号  ”  或单引号  ‘  来将变量内容结合起来,但须要特别留意, 双引号内的特殊字符可以保有变量特性,但是单引号内的特殊字符则仅为一般字符;
  5. 必要时需要以转意字符   \   来将特殊符号 ( 如 Enter, $, \, 空格符, ‘ 等 ) 变成一般符号;
  6. 在一串指令中,还需要借由 其它的shell指令 提供的信息,可以使用   ` command`  (特别特别注意,那个 ` 是键盘上方的数字键 1 左边那个按键,而不是单引号!) 或则$(命令) 。
  7. 若该变量需要扩增变量内容时,则需以双引号及 $变量名称 如: “$PATH”:/home  继续累加内容;
  8. 若该变量需要在其它子程序执行,则需要以  export  来使变量变成环境变量,例如:

#export  PATH

  1. 通常大写字符为系统预设变量,自行设定变量可以使用小写字符,方便判断 ( 纯粹依照使用者兴趣与嗜好 ) 。
  2. 取消变量的方法为: unset  变量名称。例如:

#unset   myname

环境变量

在所有 Unix类Unix系统中, 每个进程都有其各自的环境变量设置。 缺省情况下, 当一个进程被创建时, 除了创建过程中的明确更改外,它继承了其父进程的绝大部分环境设置。 在API层级上, 使用fork和exec函数进行变量设置。或利用bashshell文件, 使用特殊的命令调用来改变环境变量:通过env 间接替代或者使用ENVIRONMENT_VARIABLE=VALUE <command> 标识. 所有的Unix操作系统 以及DOSMicrosoft Windows 都是用环境变量,但是它们使用不同的环境变量名称。我们可以通过运行程序来访问环境变量的值。环境变量的例子包括:

Shell 脚本 和 批处理文件 使用环境变量来存储临时值,用于以后在脚本中引用,也用于传递数据和参数给子进程。 在Unix系统中,一个在脚本或程序中更改的环境变量值只会影响该进程,亦可能影响其子进程。其父进程和无关进程将不受影响。在DOS中,更改或删除一个批处理文件中的环境变量值将改变变量的期限命令的存在。

在Unix系统通过初始化脚本启动时,环境变量通常会在此时被初始化,因此会被系统中的其它进程所继承。用户可以而且经常添加环境变量到他们使用的shell脚本中。 在Windows系统中,环境变量的缺省值存储在 windows 注册表 中,或者在 autoexec.bat 自动执行的批处理文件中设置。(来自维基百科)

Linux是一个多用户多任务的操作系统,可以在Linux中为不同的用户设置不同的运行环境,具体做法是设置不同用户的环境变量。

Linux环境变量分类

一、按照生命周期来分,Linux环境变量可以分为两类:
1、永久的:需要用户修改相关的配置文件,变量永久生效。
2、临时的:用户利用export命令,在当前终端下声明环境变量,关闭Shell终端失效。

二、按照作用域来分,Linux环境变量可以分为:
1、系统环境变量:系统环境变量对该系统中所有用户都有效。
2、用户环境变量:顾名思义,这种类型的环境变量只对特定的用户有效。

Linux设置环境变量的方法

一、在/etc/profile文件中添加变量 对所有用户生效(永久的)
用vim在文件/etc/profile文件中增加变量,该变量将会对Linux下所有用户有效,并且是“永久的”。
例如:编辑/etc/profile文件,添加CLASSPATH变量

vim /etc/profile

export CLASSPATH=./JAVA_HOME/lib;$JAVA_HOME/jre/lib

注:修改文件后要想马上生效还要运行source /etc/profile不然只能在下次重进此用户时生效。
二、在用户目录下的.bash_profile文件中增加变量 【对单一用户生效(永久的)】
用vim ~/.bash_profile文件中增加变量,改变量仅会对当前用户有效,并且是“永久的”。

vim ~/.bash.profile

export CLASSPATH=./JAVA_HOME/lib;$JAVA_HOME/jre/lib

注:修改文件后要想马上生效还要运行$ source ~/.bash_profile不然只能在下次重进此用户时生效。
三、直接运行export命令定义变量 【只对当前shell(BASH)有效(临时的)】
在shell的命令行下直接使用export 变量名=变量值
定义变量,该变量只在当前的shell(BASH)或其子shell(BASH)下是有效的,shell关闭了,变量也就失效了,再打开新shell时就没有这个变量,需要使用的话还需要重新定义。

(原文:http://www.jianshu.com/p/ac2bc0ad3d74

变量的有效范围

变量也有使用的“范围”?在export命令说明中,就提到了这个概念。如果在跑程序的时候,有父进程与子进程的不同程序关系时,则“变量”可否被引用与export有关。被export后的变量,我们可以称它为“环境变量”。环境变量可以被子进程所引用,但是其他的自定义变量内容就不会存在于子进程中。

  • 在某些不同的书籍中会谈到“全局变量”(global variable)与“局部变量”(local variable)。基本上可以这样区分:环境变量=全局变量   自定义变量=局部变量

为什么环境变量的数据可以被子进程所引用?这是因为内存配置的关系。理论上是这样的:

  • 当启动一个shell,操作系统会分配一记忆块给shell使用,此内存内的变量可让子进程取用;
  • 若在父进程利用export功能,可以让自定义变量的内容写到上述的记忆块当中(环境变量);
  • 当加载另一个shell时(即启动子进程,而离开原本的父进程了),子shell可以将父shell的环境变量所在的记忆块导入自己的环境变量块当中。

通过这样的关系,我们就可以让某些变量在相关的进程之间存在,以帮助自己更方便地操作环境。不过要提醒的是,这个“环境变量”与“bash的操作环境”意思不太一样,举例来说,PS1并不是环境变量,但是这个PS1会影响到bash的接口(提示符)。

vi与vim

讲真,在这之前我并不知道vi与vim是有区别的,一直以为是因为偷懒将vim缩写成vi,然而并不是。

    vim是从 vi 发展出来的一个文本编辑器 。代码补完、编译及错误跳转等方便编程的功能特别丰富,在程序员中被广泛使用。和Emacs 并列成为类Unix系统 用户最喜欢的编辑器。vim是vi的加强版,它不仅兼容 vi 的所有指令,而且还有一些新的特性在里面。新的特性主要体现在以下几个方面: 1、多级撤消 2、易用性 3、语法加亮 4、可视化操作 5、对vi的完全兼容,所以在某些情况下可以把vim当做vi来使用。当然,作为生信工作者使用最多的还是vim,不说其他,语法加亮这一特性就给我们的编程带来很多便利。要想完全掌握vim的操作指令很耗脑容量,也没必要,记住一些常用指令就行。

vim指令大全:

vim的环境设置与记录(我的Vim之路