myrna本地化的安装与使用

一. Myrna的简介

Myrna主要用途是,使用RNA-seq数据来计算基因差异表达。
Myrna的运行步骤是:

1. 将reads进行预处理; 
2. 使用Bowtie将reads比对到参考序列上; 
3. 使用R/Bioconductor将比对结果定位到genes上; 
4. 标准化基因的比对数; 
5. 计算差异表达的统计数据; 
6. 将结果绘制出图。

Myrna有三种使用方式:

1. 在云端使用,Amazon's Elastic MapReduce service; 
2. 在 hadoop 集群上使用; 
3. 在单独的一台计算机上使用。

参考文献(2010):Cloud-scale RNA-sequencing differential expression analysis with Myrna
只有读该文献才能全面了解Myrna.

二. Myrna的安装

1. 下载Myrna并安装

cd /home/chenlianfu/programs/
$ wget http://nchc.dl.sourceforge.net/project/bowtie-bio/myrna/1.2.0/myrna-1.2.0.zip
$ unzip myrna-1.2.0.zip

设置Myrna的环境变量MYRNA_HOME:
$ vim ~/.bashrc
export MYRNA_HOME=/home/chenlianfu/programs/myrna-1.2.0

2. 安装bowtie

安装完bowtie,再将其执行程序bin目录添加到环境变量MYRNA_BOWTIE_HOME。或者在运行程序时候使用参数 –bowtie

设置Myrna的环境变量MYRNA_BOWTIE_HOME:
$ vim ~/.bashrc
export MYRNA_BOWTIE_HOME=/home/chenlianfu/programs/bowtie-1.0.0/

3. 安装R,Bioconductor及必须的包

本软件(版本1.20)的使用必须使用R2.14,高版本的R还不兼容。故使用Myrna自带的一个脚本来下载R2.14并安装。

# R

安装R包 lmtest 和 multicore:
> install.packages("lmtest")
> install.packages("multicore")

安装Bioconductor:
> source("http://bioconductor.org/biocLite.R")
> biocLite()

安装Bioconductor包 IRanges,geneplotter,BiocGenerics 和 biomaRt:
> source("http://bioconductor.org/biocLite.R")
> biocLite("IRanges")
> biocLite("geneplotter")
> biocLite("BiocGenerics")
> biocLite("biomaRt")

设置Myrna的环境变量MYRNA_RHOME(包含bin/R和bin/Rscript):
vim ~/.bashrc
export MYRNA_RHOME=/usr/local/

4. 安装SRA toolkit

如果输入Myrna的软件是.sra格式,则需要SRA toolkit中的程序 fastq-dump 了。安装玩 SRA toolkit 后,再设置Myrna的环境变量:

$ vim ~/.bashrc
export MYRNA_SRATOOLKIT_HOME=/home/chenlianfu/programs/sratoolkit.2.3.1-centos_linux64/bin/

5. 测试本地化运行程序

使用本地化运行程序myrna_local来进行测试,通过表示安装成功。

$ $MYRNA_HOME/myrna_local --test
...
PASSED install test

6. 设置的一些环境变量

$ vim ~/.bashrc
export MYRNA_BOWTIE_HOME=/home/chenlianfu/programs/bowtie-1.0.0/
export MYRNA_HOME=/home/chenlianfu/programs/myrna-1.2.0
export MYRNA_RHOME=/usr/local/
export MYRNA_SRATOOLKIT_HOME=/home/chenlianfu/programs/sratoolkit.2.3.1-centos_linux64/bin/
export MYRNA_REFS=/home/chenlianfu/programs/myrna-1.2.0/reftools/

三. Myrna的使用

1. 本地化运行程序的使用参数:

--reference <path>
    参考序列的文件夹。该文件夹由一个参考序列jar文件使用unzip解压缩出来的,里面主
要包括两个子文件夹:1."index":里面包含bowtie-bulid产生的index文件;2."ival
s":里面包含gene的注释信息。

--input <path>
    输入的文件或文件夹路径。如果不指定--preprocess 参数,则该参数后接的是一个
文件夹,该文件夹包含了reads的预处理结果,即 --preprocess-output 的输出结果;
如果指定了 --preprocess 参数,则是一个文本文件,内容为输入的序列文件的一个清单
(Labeled manifest files)。该文件描述了输入序列文件和样本的信息。序列文件是
FASTQ或sra格式文件,FASTQ格式可以是gzip或bzip2压缩。该文件每一行代表一个测序
的结果,是以下是 unpaired reads 和 paired reads 的书写方法:
<URL>(tab)<Optional MD5>(tab)<Sample label>
<URL 1>(tab)<Optional MD5 1>(tab)<URL 2>(tab)<Optional MD5 2>(tab)<Sample label>

其中URL是测序结果文件的路径,可以是ftp站点上的文件;MD5是可选的,以利于下载文件后
用于文件的完整度检验,若不想进行检验,则设置为0;<Sample label>的写法为:
<Group ID>-<BioRep ID>-<TechRep ID>

其中,Group ID 是样品的名称;BioRep ID是生物学重复的代号;TechRep ID是技术性
重复的代号。比如:
Tumor-Subject1-Lane1
Tumor-Subject1-Lane2
Tumor-Subject2-Lane1
Tumor-Subject2-Lane2
Normal-Subject1-Lane1
Normal-Subject1-Lane2
Normal-Subject2-Lane1
Normal-Subject2-Lane2

--output <path>
    输出文件的路径。如果 --just-preprocess 参数给定,则只是对reads文件进行预
处理,给出结果。默认情况下,则是输出最终结果:一系列的表格,图和比对文件。

--intermediate <path>  default: /tmp/myrna/intermediate/<PID>
    中间结果文件存放的临时路径。如果指定了 --keep-intermediates 或 --keep-
all参数,则该文件则会一直存在。默认情况下会使用完毕后删除。

--preprocess-output <path>
    对reads进行预处理的输出路径。由于对reads的预处理会消耗很多时间,因此保留数据
的预处理结果,对再次运行程序来说,可以节约很多时间。

--keep-intermediates <path>
    保留程序运行中生成的结果文件夹和文件,默认情况下这些文件会被尽快的自动删除。

--keep-all
    保留所有的临时文件,默认情况下这些文件会被尽快的自动删除。

--cpus <int>  default: 1
   设定程序运行使用的CPU线程数。

--bowtie <path>
    设定Bowtie的路径,以运用bowtie进行比对的步骤,该参数优先于环境变量MYRNA_
BOWTIE_HOME,$MYRNA_HOME/bin的子目录,以及 PATH。

--fastq-dump <path>
    设定 SRA toolkit 中的 fastq-dump 程序的位置,优先于PAHT。

--Rhome <path>
    设定R和R/Bioconductor安装程序的位置,该目录包括 bin/R 和 bin/Rscript,
优先于环境变量MYRNA_RHOME和PATH。

2. emr, hadhoop 和 local运行的共同参数

--quality { phred33 | phred64 | solexa64 }  default: phred33
    设置输入reads的碱基质量格式。phred33:输入的碱基质量等于ASCII码值加上33;最
低碱基质量是“#”. phred64:输入的碱基质量等于ASCII码值加上64;最低碱基质量是“B”. 

--preprocess
    当 --input 的参数是一个文本文件(含有输入序列信息的清单)时,则必须加入该参数
;该参数指定出了文本文件中的序列信息来,Myrna则进行这些reads的预处理,将生成的预处
理reads结果输出到 --preprocess-output 指定的文件夹;如果不指定--preprocess
-output参数,则输出到 --intermediate 指定的文件夹。

--just-preprocess
    和 --preprocess 参数类似,但不同的是进行reads预处理后,即停止pipeline的
运行,并将结果输出到 --output指定的文件夹。

--just-align
    不将Myrna的pipeline运行到底,只是运行到比对结束,并将结果存放到 --output 
参数指定的URL.

--resume-align
    使用此命令来运行上一个命令(比对完毕)过后的阶段。 此时 --input 的参数则是上
一个命令中 --output的URL.

--bowtie-args "<args>"  default: -m 1
    设置bowtie的参数。

--discard-reads <fraction>  default: 0.0
    丢弃一定比例输的reads来进行程序的运行。比如 0.5 则代表丢弃50%的reads。该参
数应用于所有输入的reads。对于程序的调试比较有用。

--influence <int>  default: 1
    一个read的"influence"长度。和下面3个参数有关系。

--from3prime
    设定修剪过后的reads从3'端开始衡量其"influence"。此参数是默认设置,和下2个
参数是冲突的,三选一。

--from5prime
    设定修剪过后的reads从5'端开始衡量其"influence"。

--from-middle
    设定修剪过后的reads从中间开始衡量其"influence"。

--top <int>  default:50
    将genes按p值从小到大排序,Myrna将报告出前一定数目的genes的比对结果。默认是
50个genes。

--family { poisson | gaussian }  default: poisson
    设定检测gene的差异表达的统计方法。possion适合于样本<=10个;gaussian适合于
样本数>10。当然这不是硬性规定,具体还要用实验验证。当 --nulls 参数是 0 或者 没有
给出时,则使用此参数中的模型来作参数检验(parameric statistical model),计算
出p-value.当 --nulls参数大于0,则使用permutation test来计算p-value.

--normalize { total | median | upper-quartile }
    Each count (or pooled count) is "normalized" according to a 
per-label baseline normalization factor.

--gene-footprint { union | intersect }  default: intersect
    计算比对和genges的关联时所使用的 gene footprint。intersect 和 union 
分别对应这reference文件夹中的ivals子文件夹中的ui和un文件夹;其中un对应的基因区
要比ui多,并且un的基因区包含了ui的,个人推测可能是un包含了5'端和3'端非翻译区。
union: 将一个比对结果指向某一个gene的必要条件是,read的"influence"(默认下是
从3'端开始)有一个碱基和一个gene任意的(any)transcripts重叠;intersect: 将
一个比对结果指向某一个gene的必要条件是,read的"influence"(默认下是从3'端)有一
个碱基和一个gene所有的(all)的transcripts重叠。我个人的理解,通俗来讲:以上数个
参数的目的是计算比对到gene上的reads的数目。根据比对结果,默认将reads从3'端开始
计算,如果比对上reads有一个碱基和一个gene的一个外显子重叠,则在union模式下这个
gene则可以得到一个read比对结果;而在intersect模式下,这个reads必须和该gene所
有的外显子有一个碱基的重叠,则该gene才能得到一个read比对结果。可见intersect要求
read贯穿整个gene才行能得到比对到gene的数目;因此,对于short reads进行基因表达
量计算的时候,则需要选union模式,不能选默认的模式。

--paired-ttest
    如果设定此值,则表示: 1, 是对2个样作比较; 2, 这2个样的生物学重复必须要匹
配。当该参数使用时,将对两个样作t检验。

--nulls <int>
    Set the number of null statistics generated per gene to <int>.
 A null statistic for a gene is calculated by permuting the sample
 labels for the counts and normalization factors, then re-calcula-
-ting the statistical test. 如果<int>是一个非0的数,则p-values由permu-
-tation test得出;否则,p-values由参数检验模型(parametric statistical 
model)计算出。

--truncate <int>
    将所有的reads都从3'端开始修剪,直到该参数设定的长度;小于该长度的reads依然
保留。

--truncate-discard <int>
    将所有的reads都从3'端开始修剪,直到该参数设定的长度;而小于该长度的reads全
部剔除。

--ditch-alignments
    不输出比对结果和图形结果;依然报告差异表达基因,p-values和相关的图形结果。

--discard-mate <int>
    对于 paired-end reads,去掉其中一端的序列( 1 或 2 )再来进行比对。这对于
一些数据集的比较有用,比如数据集中同时有 paired reads 和 unpaired reads时。

--pool-reps
    将所有的重复,包括生物学重复和技术性重复融合后,再进行计算统计。

--pool-tech-reps
    将技术性重复进行融合后,在进行计算统计。

--test
    搜索Myrna需要的工具(Bowtie, R/Bioconductor 和 fastq-dump),来检测
Myrna能否正确运行。

--tmpdir `<int>`  default: /tmp/Myrna/invoke.scripts
    临时文件存放的位置,是一个动态生成的脚本。

四. 个人感受

使用时,生成的结果没有rpkm值,p-value值及预期的图表结果。不知到哪儿出错了。比对的之前的过程都是正常的。同时使用sra文件时候,提示fastq-dump找不到或无法执行(实际上在–test的时候可行)。无法正常使用myrna。费心费力,不想深究了…悲剧…再去研究cufflinks去。

使用代理上国外网站

经常的,会遇到重要的国外网站上不去。比如著名的bowtie软件,在国内就上不去该网站了。或者google时常抽风。欲解决这些问题,则要开启一些代理来正常访问国外的网站,以下是具体的方法:

1. 从一些网站搜索大量的代理IP

代理IP网来得到一些代理IP。此时打开该网站主页,选择国外代理页面。

$ wget http://www.dlipw.com/a/guowaidailiip/2013/0330/4269.html
下载了该网页的所有内容
$ grep -P "\d+\.\d+\.\d+\.\d+:\d+" 4269.html > IPs
通过perl的正则表达式提取出含ip地址的行
$ perl -p -i -e 's/^.*(\d+\.\d+\.\d+\.\d+:\d+).*$/$1/' IPs
继续将ip地址和端口提取出来,每行一个ip地址和端口信息

2. 代理猎手软件Proxy Hunter来验证IP

代理猎手软件解压缩后,在windows下打开直接使用。代理猎手下载网页:http://www.orsoon.com/soft/down3555.html。其使用方法:

验证数据设置:
选择系统菜单——参数设置——验证数据设置——添加——验证地址输入:http://bowtie-bio.s
ourceforge.net/index.shtml——特征字串输入:ultrafast——确定。
导入并验证IP:
点击搜索结果——导入结果,通过浏览来打开上一步骤中的文件——验证全部——Free的结果即为我
们需要的代理IP了。

使用以上方法搜索出的几个结果:

2.133.93.50:9090
2.135.237.92:9090
2.135.237.98:9090
2.181.177.7:8080
2.183.155.2:8082
2.184.30.2:8080
2.184.31.2:8080
5.8.242.11:8080
2.184.30.2:8080
2.183.155.2:8082
5.199.166.250:3128
1.34.181.22:110

3. firefox方便地代理上网

安装附加组件AutoProxy方便地代理上网,以免自己手动切换的麻烦。该组件设置一些规则,指定一些网站代理上网,其它网张正常上网。

4. 一些在线网页版的代理

http://proxyie.cn/

mplayer播放参数与设置

一. mplayer双字幕设置

直接将中英字幕文件合并成一个文件即可,即可直接用于双字幕mplayer播放。

$ cat my_film.chs.srt my_film.eng.srt > my_film.srt
$ mplayer -ass my_film.mkv -sub my_film.srt
当中文字母加载出现乱码,使用 -ass 参数。

二. mplayer播放调节参数

1. 字幕参数:

x 和 z
    调整字幕延迟增加/减少 0.1 秒
r 和 t
    上/下调整字幕位置
v
    关闭或打开字幕
j
    切换字幕,对应的字幕位于视频同一目录下并文件名一致

2. 声音参数:

+ 和 -
    调整音频延迟增加/减少 0.1 秒
/ 和 *
    降低/提高音量
9 和 0
    降低/提高音量
m
    静音开关
#
    循环可用的音轨

3. 播放控制参数:

左方向键 和 右方向键
    后退/快进 10 秒钟
上方向键 和 下方向键
    快进/后退 1 分钟
pgup 和 pgdown
    快进/后退 10 分钟
[ 和 ]
    减少/增加当前播放速度 10%
{ 和 }
    减半/加倍当前播放速度
p / SPACE
    暂停 (再按取消暂停)
q / ESC
    停止播放和退出

4. 屏幕参数

f
    切换全屏
T
    切换置顶
o
    循环OSD状态:即显示播放计时
P
    在 OSD 上显示进度条、已播放时间以及总长度信息

PubMed与PMC的区别

PubMed Central (PMC 公共医学中心)是一个提供有关生命科学与生物医学的回溯性电子期刊全文数据库,它是在2000年1月由隶属美国国立图书馆(NLM)的国家生物技术信息中心(NCBI)所创建与管理的。PMC采取自愿加入的原则,某期刊一旦加入,必须承诺期刊出版后一定时期内(最好六个月,不超过1年)将其全文提交给PMC,由PMC提供免费全文检索和访问。目前已收录5百种期刊。 PMC与PubMed的关系:两者都是NLM建立的数据库。其中PubMed是一个基于互联网的文献检索系统,它收录了几千种生命科学期刊的目次和文摘,目前数据已回溯至1966年,该数据库提供了与PMC全文的链接以及与数千种期刊网站的链接。而PMC是由NLM建立的免费生命科学电子期刊全文数据库,目前收录期刊百余种,PMC的所有论文在PubMed中都有相应的记录。

Quake的安装和使用

一. Quake简介

Quake是由CBCB(Center for Bioinformatics and Computational Biology)开发的运用于修正序列错误的软件。Quake采用k-mer的错误修正方式,特别适合于Illumina测序的short reads数据,将reads中的错误碱基进行修正,同时,必须要满足碱基的基因组覆盖度要>15X。其文章2010年发表在Genome Biology上。

二. Quake的安装

1. 下载并安装Boost
2. 下载Quake并安装

$ wget http://www.cbcb.umd.edu/software/quake/downloads/quake-0.3.4.tar.gz
$ tar zxf quake-0.3.4.tar.gz
$ cd Quake/src
$ make

3. 安装JELLYISH,并将jellyfish链接到quake

$ ln -s ..../jellyfish ../bin/jellyfish
这一步需要安装jellyfish,然后将其软链接到Quake的bin目录下。或者修改Quake的bin
目录下的quake.py脚本,将jellyfish所在的目录进行修正。

4. 安装R及R的软件包VGAM

$ R
> install.packages("VGAM")
> q(save="no")

5. 使用Quake的bin目录下的quake.py来运行程序

三. quake的使用参数

1. 主要参数:

-r READSF
    Fastq文件名
-f READS_LISTF
    一个文件,其中包含fastq文件名,每行一个文件名;如果是双末端reads,则是一行两
个文件名
-k K
    使用的k-mer的长度。如果基因组大小为G,则k-mer长度选择为: k ~= log(200G)
/log(4)
-p PROC  default: 4
    使用的CPU线程数
-q QUALITY_SCALE
    使用的碱基质量格式,一般是64或33.如果不给出,则软件会自行猜测

2. 计算k-mer参数:

--no_jelly
    使用quake自带的一个简单的程序来进行k-mer计数,而不是使用jellyfish。该程序
是单线程运行,速度相对慢,能满足像微生物基因组这样的小基因组测序,但是不适合于大的基
因组。
--no_count  default: false
    不进行k-mer计数,而其计数结果已经存于在目标文件[readsfile].qcts和[reads
file].cts中了。
--int  default:false
    对kmers以整数的方式进行计数,而不使用碱基质量值。默认情况下是利用到了碱基质量,
计数的称为qmer。
--hash_size=HASH_SIZE
    Jellyfish的参数,用来设置Hash的大小。如果不设置,Quake则会使用k来估计出一
个值来。

3. 覆盖度模型参数:

--no_cut  default: false
    Quake使用k-mer计数来画直方图,通过调用R脚本cov_model_qmer.r来画直方图,并
决定出Coverage cutoff的阈值。默认情况下是最优化的模型,从而得出一个cutoff的Co
verage值,此值将输出到文件cutoff.txt中。
--ratio=RATIO
    确定Coverage cutoff值的时候,该值处对应的qmer错误的可能性是正确的可能性的
倍数(默认是200倍),即正确率不足0.5%。该值越小,则阈值越松。

4. reads的修正参数:

-l MIN_READ
    输出的长度>=此值的修正后的reads。很多reads由于修正或修剪后会较短
--headers
    输出的fastq文件使用的头和原始文件一致,不包含修正信息
--log
    输出一个log文件,里面包含所有的修正日志,包括"碱基质量 位置 新的碱基 旧的碱基"。

四. Quake的结果

使用Quake对双末端测序文件 A.fastq 和 B.fastq 进行reads修正,则产生如下结果文件:

A.cor.fastq  B.cor.fastq
    修正过后的reads文件               
A.cor_single.fastq  B.err_single.fastq
    A文件中修正过后的reads结果;而B文件中的序列则是error reads
A.err_single.fastq  A.cor_single.fastq
    B文件中修正过后的reads结果;而A文件中的序列则是error reads
A.err.fastq  B.err.fastq
    A和B文件中成对的reads都是error reads

五. Quake的correct程序

Quake的bin文件中有一支C++的程序,用于reads的修正。其实Quake.py在运行过程会调用该程序,但是当出现基因组覆盖度较低,而cutoff的值出现异常结果,比如小于1时,可以考虑使用该程序来重新对reads进行修正,取cutoff值为1.0。其使用方法为:

$ correct -f [fastq list file] -k [k-mer size] -c [cutoff] -m [counts file] -p 24

jellyfish的安装和使用

一. JELLYFISH简介

JELLYFISHCBCB(Center for Bioinformatics and Computational Biology)的Guillaume Marçais 和 Carl Kingsford 研发的一款计数 DNA 的 k-mers 的软件。该软件运用 Hash 表来存储数据,同时能多线程运行,速度快,内存消耗小。该软件只能运行在64位的Linux系统下。其文章于2011年发表在杂志 Bioinformatics 上。

二. JELLYFISH安装

$ wget http://www.cbcb.umd.edu/software/jellyfish/jellyfish-1.1.10.tar.gz
$ tar zxvf jellyfish-1.1.10.tar.gz
$ mkdir jellyfish
$ cd jellyfish-1.1.10
$ ./configure --prefix=$PWD/../jellyfish
如果安装在当前目录中,会报错。
$ make -j 8
$ make install

三. JELLYFISH的使用

1. jellyfish的使用方法

jellyfish的功能有:kmer计数;融合二进制的Hash结果;统计Hash结果;通过Hash结果来画直方图;将Hash结果输出成文本格式;查询指定k-mer的数目。

$ jellyfish count [-o prefix] [-m merlength] [-t threads] [-s hashsize] [--both-strands] fasta [fasta ...]
$ jellyfish merge hash1 hash2 ...
$ jellyfish dump hash
$ jellyfish stats hash
$ jellyfish histo [-h high] [-l low] [-i increment] hash
$ jellyfish query hash
$ jellyfish cite

2. k-mer的计数

使用 count 的命令来执行计数功能,使用例子:

$ jellyfish count -m 16 -s 100M -t 24 -o mer_counts -c 7 input.fastq
使用fastq文件在默认参数上和fasta文件没有区别。生成的hash结果为二进制文件。

常用参数:

-m | --mer-len=<num>
    使用的k-mer的长度。如果基因组大小为G,则k-mer长度选择为: k ~= log(200G)
/log(4)。
-s | --size=<num>
    Hash 的大小。最好设置的值大于总的独特的(distinct)k-mer数,这样生成的文件只
有一个。若该值不够大,则会生成多个hash文件,以数字区分文件名。如果基因组大小为G,每
个reads有一个错误,总共有n条reads,则该值可以设置为『(G + n)/0.8』。该值识别 
M 和 G。
-t | --threads=<num>  default: 1
    使用的CPU线程数
-o | --output=<string>  default: mer_counts
    输出的结果文件前缀
-c | --counter-len=<num>  default:7
    k-mer的计数结果所占的比特数,默认支持的最大数字是2^7=128。对于基因组测序覆盖
度为N,则要使设置的该值要大于N。该值越大,消耗内存越大。
-out-counter-len=<num>  default:4
    输出的二进制hash文件中的计数结果所占的字节数,一个字节是8比特。则默认支持的最大
数字是2^32=4.3G
-C | --both-strand  default: false
    对正义链和反义链都进行计数
-q | --quake  default: false
    quake兼容模式
--quality-start=<num>  default: 64
    起始碱基质量的ASCII值,默认为PHRED64
--min-quality=<num>  default: 0
    支持的最小的碱基质量值,低于此值的碱基将由N代替
-L | --lower-count=<num>
    不输出数目低于此值的k-mer
-U | --upper-count=<num>
    不输出数目高于此值的k-mer

3. 融合二进制的输出结果

上一步的输出结果为二进制文件,可能输出了多个hash文件,因此需要将这些hash文件合并成一个文件,此时用到 merge 命令。使用方法:

$ jellyfish merge -o mer_counts_merged.jf hash1 hash2 ...

常用参数:

-o | --output=<string>  default: mer_counts_merged.jf
    输出的结果文件
--out-counter-len=<num>  default: 4
输出的二进制hash文件中的计数结果所占的字节数,一个字节是8比特。则默认支持的最大数字
是2^32=4.3G

4. 对hash结果进行统计

k-mer的结果以hash的二进制文件结果给出,需要统计出k-mer总数,特异的k-mer数目,只出现过一次的kmer数,出现了最多的k-mer的数目等信息。以stats命令来运行。使用方法:

$ jellyfish stats hash
示例结果为:
Unique:    32355544    #只出现过一次的k-mer的数目
Distinct:  88414020    #特异性的k-mer数目,包含上一个的数据
Total:     432232807   #总的k-mer数目
Max_count: 85348       #同一个k-mer出现的最多的数目

常用参数:

-L | --lower-count=<num>
    不统计数目低于此值的k-mer
-U | --upper-count=<num>
    不统计数目高于此值的k-mer

5. 通过Hash结果来画直方图

对k-mer的计数结果有个直观的认识,则需要统计出现了x(x=1,2,3…)次的kmer的数目y,以x,y为横纵坐标画出直方图。使用 histo 命令能给出 x 和 y 对应的值,将结果默认输出到标准输出。其使用方法为

$ jellyfish histo -l 1 -h 1000 hash

常用参数:

-l | --low=<num>  default: 1
    最低的 x 轴的值。同时结果会将低于此值的所有的k-mer的数目作为 (x-1) 的值。因
此该值为 2 和 1 的结果是一致的。
-h | --high=<num>  default: 10000
    最高的 x 轴的值。同时结果会将高于此值的所有的k-mer的数目的和作为 (x+1) 的值。
-i | --increment=<num>  default: 1
    x 轴取值是每隔该数值取值
-t | --threads=<num>  default: 1
    使用的CPU线程数
-f | --full  default: false
    全部的直方图

6. 将二进制Hash结果转换成文本文件

由于count命令生成的结果为二进制的,如有需要,则可以转换成可读文本文件。使用 dump 命令,使用方法:

$ jellyfish dump -c -t -U 1000 hash

常用参数:

-c | --colum  default: false
    生成结果为2列,第一列为k-mer序列,第二列为对应的数目。默认情况下是是fasta格
式,fasta的头为k-mer的数目,fasta的序列为k-mer的序列。
-t | --tab  default: false
    当 -c 参数存在时,以tab来进行分隔两行。默认是以空格来分开的。
-L | --lower-count=<num>
    不输出小于该值的k-mer
-U | --upper-count=<num>
    不输出高于该值的k-mer
-o | --output=<file>
    输出文件的路径和名称

7. 查询指定的k-mer出现的次数

如果需要从Hash结果中查询指定的k-mer出现的次数,则要是用 query 命令。从标准输入读取k-mer的序列,从标准输出得到k-mer对应的数目。使用方法

$ jellyfish query hash

常用参数:

-C | --both-strands  default: false
    同时查询k-mer序列的正负链
-i | --input=<file>
    输入的文件
-o | --output=<file>
    输出的文件

四. 思考

1. 对Illumina paired-end测序结果进行jellyfish分析

由于paired-end序列有一定的顺序,需要将第2个文件的序列进行反向重复后,在和第一个文件的序列合到一起进行分析。可以使用Trinity中附带的软件fastool来将fastq文件转换成fasta文件,以及反向重复的转换。

$ $Trinity_Home/trinity-plugins/fastool/fastool --illumina-trinity --to-fasta reads_1.fastaq > reads_1.fasta
$ $Trinity_Home/trinity-plugins/fastool/fastool --rev --illumina-trinity --to-fasta reads_2.fastaq > reads_2.fasta
$ cat reads_1.fasta reads_2.fasta > both.fasta
$ jellyfish count ....

mail的使用

1. mail的使用使用的是mail命令,或mailx命令。

使用 which mail 查看到mail命令的所在,发现其实是指向 mailx 命令。
使用参数:

-s subject
    邮件的标题。若标题有空格,则需要使用引号。
-a attachment
    将目标文件作为附件发送。若有多个附件需要发送,则使用多个该参数。
-c address
    抄送副本到邮件地址列表。这些邮件地址使用逗号分隔。抄送的邮件地址和收件人地址能
被所收件地址看到。
-b address
    暗送的邮件地址列表。这些邮件地址使用逗号隔开。暗送的邮件地址不能被其收件地址看
到。故mail命令不能将邮件分别发送到邮件地址列表。

使用方法:

$ mail -s "a e-mail subject" -a ./test.tar.gz chenllianfu@foxmail.com < mail_content
$ cat mail_content | mail -s "a e-mail subject" -a ./test.tar.gz chenllianfu@foxmail.com
$ echo "mail_content" | mail -s "a e-mail subject" -a ./test.tar.gz chenllianfu@foxmail.com
$ mail -s "a e-mail subject" -a ./test.tar.gz chenllianfu@foxmail.com
input
EOT

2. 提高默认允许的附件大小

If you get the message:

postdrop: warning: uid=500: File too large
sendmail: fatal: : message file too big

You can fix this bij editing /etc/postfix/main.cf and add:

message_size_limit = 52428800

(example creates a limit of 50Mb)

tar命令的使用

tar命令主要用来对数据进行压缩和打包。其常用的参数有:

-C | --directory=DIR
    改变目录至DIR。在打包过程中,不指定此参数,则压缩包可能会包含很长很长的路径,解
压后会生成一长串的文件路径。
--remove-files
    在添加文件至归档后,然后再删除这些备份过了的文件,以节约磁盘空间。
-f | --file=ARCHIVE
    生成的归档文件。将生成的归档文件指定到具体的路径和名称,与 -C 参数无关。
-c | --create
    创建一个新的归档。
-x | --extract
    从归档中解出文件。
-t | --list
    列出归档卷标。
-z | --gzip,--ungzip
    通过gzip来过滤归档
-j | --bzip2
    通过bzip2来过滤归档
-J | --xz
    通过xz过滤归档
-Z | --compress, --uncompress
    通过compress过滤归档