使用HTSeq进行有参转录组的表达量计算

1. HTSeq简介

HTSeq是使用Python编写的一支用于进行基因Count表达量分析的软件,能根据SAM/BAM比对结果文件和基因结构注释GTF文件得到基因水平的Counts表达量。HTSeq进行Counts计算的原理非常简单易懂,容易上手。

2. HTSeq安装

PYPI下载HTSeq的Python包
$ wget https://pypi.python.org/packages/46/f7/6105848893b1d280692eac4f4f3c08ed7f424cec636aeda66b50bbcf217e/HTSeq-0.7.2.tar.gz
$ tar zxf HTSeq-0.7.2.tar.gz
$ cd HTSeq-0.7.2
$ /opt/sysoft/Python-2.7.11/bin/python setup.py build
$ /opt/sysoft/Python-2.7.11/bin/python setup.py install
$ cd ../ && rm HTSeq-0.7.2 -rf

3. HTSeq使用

3.1 HTSeq的Count模式

HTSeq计算counts数有3种模式,如下图所示(ambiguous表示该read比对到多个gene上;no_feature表示read没有比对到gene上):
HTSeq Count模式

3.2 HTSeq的使用命令

HTseq安装完毕后,在Python软件的bin目录下生成htseq-count命令。
htseq-count运行简单示例:

对于非链特异性真核转录组测序数据
$ /opt/sysoft/Python-2.7.11/bin/htseq-count -f sam -r name -s no -a 10 -t exon -i gene_id -m union hisat2.sam genome.gtf > counts_out.txt
对于链特异性测序真核转录组测序数据
$ /opt/sysoft/Python-2.7.11/bin/htseq-count -f sam -r name -s reverse -a 10 -t exon -i gene_id -m union hisat2.sam genome.gtf > counts_out.txt
对于非链特异性原核生物转录组测序数据
$ /opt/sysoft/Python-2.7.11/bin/htseq-count -f sam -r name -s no -a 10 -t exon -i gene_id -m intersection-strict bowtie2.sam genome.gtf > counts_out.txt

htseq-count命令的常用参数:

-f | --format     default: sam
  设置输入文件的格式,该参数的值可以是sam或bam。
-r | --order     default: name
  设置sam或bam文件的排序方式,该参数的值可以是name或pos。前者表示按read名进行排序,后者表示按比对的参考基因组位置进行排序。若测序数据是双末端测序,当输入sam/bam文件是按pos方式排序的时候,两端reads的比对结果在sam/bam文件中一般不是紧邻的两行,程序会将reads对的第一个比对结果放入内存,直到读取到另一端read的比对结果。因此,选择pos可能会导致程序使用较多的内存,它也适合于未排序的sam/bam文件。而pos排序则表示程序认为双末端测序的reads比对结果在紧邻的两行上,也适合于单端测序的比对结果。很多其它表达量分析软件要求输入的sam/bam文件是按pos排序的,但HTSeq推荐使用name排序,且一般比对软件的默认输出结果也是按name进行排序的。
-s | --stranded     default: yes
  设置是否是链特异性测序。该参数的值可以是yes,no或reverse。no表示非链特异性测序;若是单端测序,yes表示read比对到了基因的正义链上;若是双末端测序,yes表示read1比对到了基因正义链上,read2比对到基因负义链上;reverse表示双末端测序情况下与yes值相反的结果。根据说明文件的理解,一般情况下双末端链特异性测序,该参数的值应该选择reverse(本人暂时没有测试该参数)。
-a | --a     default: 10
  忽略比对质量低于此值的比对结果。在0.5.4版本以前该参数默认值是0。
-t | --type     default: exon
  程序会对该指定的feature(gtf/gff文件第三列)进行表达量计算,而gtf/gff文件中其它的feature都会被忽略。
-i | --idattr     default: gene_id
  设置feature ID是由gtf/gff文件第9列那个标签决定的;若gtf/gff文件多行具有相同的feature ID,则它们来自同一个feature,程序会计算这些features的表达量之和赋给相应的feature ID。
-m | --mode     default: union
  设置表达量计算模式。该参数的值可以有union, intersection-strict and intersection-nonempty。这三种模式的选择请见上面对这3种模式的示意图。从图中可知,对于原核生物,推荐使用intersection-strict模式;对于真核生物,推荐使用union模式。
-o | --samout 
  输出一个sam文件,该sam文件的比对结果中多了一个XF标签,表示该read比对到了某个feature上。
-q | --quiet
  不输出程序运行的状态信息和警告信息。
-h | --help
  输出帮助信息。

3.3 HTSeq使用注意事项

HTSeq的使用有如下注意事项,否则得到的结果是错误的:

1. HTSeq是对有参考基因组的转录组测序数据进行表达量分析的,其输入文件必须有SAM和GTF文件。
2. 一般情况下HTSeq得到的Counts结果会用于下一步不同样品间的基因表达量差异分析,而不是一个样品内部基因的表达量比较。因此,HTSeq设置了-a参数的默认值10,来忽略掉比对到多个位置的reads信息,其结果有利于后续的差异分析。
3. 输入的GTF文件中不能包含可变剪接信息,否则HTSeq会认为每个可变剪接都是单独的基因,导致能比对到多个可变剪接转录本上的reads的计算结果是ambiguous,从而不能计算到基因的count中。即使设置-i参数的值为transcript_id,其结果一样是不准确的,只是得到transcripts的表达量。

3.4 HTSeq的结果

HTSeq将Count结果输出到标准输出,其结果示例如下:

gene00001	0
gene00002	9224
gene00003	880
...
gene12300	1043
gene12301	200
__no_feature	127060
__ambiguous	0
__too_low_aQual	4951
__not_aligned	206135
__alignment_not_unique	0

使用AWStats对网站流量进行统计

1. 安装AWStats

# yum install awstats

2. 使用AWStats

常用示例:
# /var/www/awstats/awstats.pl --config=localhost.localdomain -update -output > /var/www/awstats/index.html

--config=virtualhostname
    该参数用于导入配置文件。配置文件位于 /etc/awstats 或 /usr/local/etc/awstats 目录。程序会导入awstats.virtualhostname.conf或awstats.conf配置文件。
-update
    对数据统计结果进行更新。
-output
    输出HTML结果文件。

CentOS系统NAT共享上网

现在有服务器A通过PPPOE联网,而服务器B直接和服务器A使用网线连接。若需要使B能正常上网,则需要将A的网络共享给B。此外,服务器A和B都具有多网口,并都是CentOS系统。将服务器A的网络共享给B,其对两台服务器的设置如下:

1. 服务器A的设置

1.1 服务器A的第一个网口进行pppoe连接

将网线插入到服务器A的第一个网口eth0,然后设置服务器A的PPPOE连接:

安装pppoe软件
# yum install rp-pppoe
配置pppoe设置,填写上网账号和密码,该pppoe配置名称为ppp0,保证对应的网口为eth0,设置。
# pppoe-setup
关闭ppp0的连接
# ifdown ppp0
开启ppp0的连接
# ifup ppp0
若发现ppp0连接不上,输入下面命令后再连接
# pppoe-stop

1.2 服务器A的第二个网口的IP设置

再将服务器A的eth1口和服务器B的eth1口进行连接。对服务器A的eth1口进行设置:

# setup
通过setup命令配置eth1的IP,设置其:
IP地址:192.168.1.1
子网掩码:255.255.255.0
网关:192.168.1.1
DNS1:211.69.143.1
DNS2:8.8.8.8

其中DNS1是我们学校提供的DNS服务器网址,DNS2是google提供的DNS网址。可能在不同的地方其DNS网址不一样。
然后,修改 eth0 和 eth1 的配置文件,设置这两个网口开机启动:
# vi /etc/sysconfig/network-scripts/ifcfg-eth0
# vi /etc/sysconfig/network-scripts/ifcfg-eth1
ONBOOT=yes

1.3 将服务器A的ppp0网络进行NAT共享

首先,修改配置文件/etc/sysctl.conf的一个参数来开启使用NAT进行IP转发。

$ perl -p -i -e 's/net.ipv4.ip_forward =.*/net.ipv4.ip_forward = 1/' /etc/sysctl.conf

再使用iptables命令来对ppp0进行网络共享。

生成文件 /usr/local/bin/ishare 并使其可执行,执行该命令,即可共享 ppp0 网络。
# echo '#!/bin/bash
## Internet connection shating script
sysctl -w net.ipv4.ip_forward=1
sysctl -p
iptables -X
iptables -F
iptables -t nat -X
iptables -t nat -F
iptables -I INPUT -m state --state RELATED,ESTABLISHED -j ACCEPT
iptables -I FORWARD -m state --state RELATED,ESTABLISHED -j ACCEPT
iptables -t nat -I POSTROUTING -o ppp0 -j MASQUERADE' > /usr/local/bin/ishare
# chmod 755 /usr/local/bin/ishare
# /usr/local/bin/ishare

若需要开机则启动该命令:
# echo '/usr/local/bin/ishare' >> /etc/rc.d/rc.local

1.4 重启服务器A的网络设置

服务器A的配置修改完毕,然后重启服务器A的网络,使配置生效:

# /etc/init.d/network restart
保证看到pppoe和eth1网络的正常启动。

2. 服务器B的设置

配置服务器B的 eth1 网口:

# setup
通过setup命令配置eth1的IP,设置其:
IP地址:192.168.1.2
子网掩码:255.255.255.0
网关:192.168.1.1
DNS1:211.69.143.1
DNS2:8.8.8.8

修改 eth1 的配置文件,设置这该网口开机启动:
# vi /etc/sysconfig/network-scripts/ifcfg-eth1
ONBOOT=yes

再重启重启服务器B的网络设置
# /etc/init.d/network restart
保证看到eth1网络的正常启动。

3. 会出现的问题

以上设置完毕后,通过 ping 命令来检测服务器B是否能联网。若服务器ping不通192.168.1.1,则表示服务器B和服务器A连接不通;若能ping通192.168.1.1而不能联外网,则表示服务器A没有开启共享或服务器A没有联网。需要按照上述教程逐步排查。

mycology_学科杂志分区

生物大学科分区	MYCOLOGY小学科分区	杂志名称	2014-2015影响因子
1区	1区	STUDIES IN MYCOLOGY	13.25
2区	2区	FUNGAL DIVERSITY	6.221
2区	2区	PERSOONIA	5.3
3区	2区	MYCORRHIZA	3.459
3区	3区	Fungal Ecology	2.929
3区	3区	FUNGAL GENETICS AND BIOLOGY	2.587
4区	3区	MYCOLOGIA	2.471
4区	4区	Fungal biology	2.342
3区	3区	MEDICAL MYCOLOGY	2.335
4区	4区	MYCOSES	2.239
3区	3区	World Mycotoxin Journal	2.157
4区	4区	MYCOLOGICAL PROGRESS	1.913
4区	4区	MYCOPATHOLOGIA	1.528
4区	4区	CRYPTOGAMIE MYCOLOGIE	1.524
4区	4区	LICHENOLOGIST	1.454
4区	4区	Mycoscience	1.418
4区	4区	Revista iberoamericana de micología : órgano de la Asociación Espa ola de Especialistas en Micología1.056
4区	4区	SYDOWIA	1.021
4区	4区	INTERNATIONAL JOURNAL OF MEDICINAL MUSHROOMS	0.927
4区	4区	MYCOTAXON	0.705
4区	4区	JOURNAL DE MYCOLOGIE MEDICALE	0.573

blast进行重复序列屏蔽

1. 构建数据库的时候屏蔽参考序列的重复

segmasker 屏蔽氨基酸的低复杂序列
dustmasker 屏蔽核算序列的低复杂序列
windowmasker 按照序列重复的次数来屏蔽
convert2blastmask 根据小写字母来屏蔽

这几个都可以先得到一个含有屏蔽信息的文件。然后进行 makeblastdb 的时候输入这个文件,就可以相应的 masked 数据库了。

参考:http://www.ncbi.nlm.nih.gov/books/NBK279681/

2. 比对的时候对query序列的重复进行屏蔽

blast 比对的时候,可以对 query 序列进行屏蔽。 这几个参数估计这样理解:
-seg blastp的参数,是否对query 序列使用 segmasker来屏蔽低复杂重复,默认 no
-dust blastn的参数,是否对query 序列使用 dustmasker来屏蔽低复杂重复,默认 no
-lcase_masking 对query序列的小写部分进行屏蔽
-soft_masking 是否进行软屏蔽。软屏蔽则是不会使用屏蔽的序列进行种子比对,但是可以延长时候比对。硬屏蔽,则是直接不对屏蔽序列部分进行比对。blastn的默认值是yes,blastp的默认值是no

文档编辑经验点

1. 分节符的使用
点击:“页面布局”——“分隔符”——“分节符下一页”,在指定位置插入分节符,用于将文章不同的章节进行分割。这样可以保证:下一章节的第一行则总是在页面的最上面;下一章节的排版和上一章节可以不一致,例如纸张方向不一致。

2. 使用Endnote分别对每一章节插入文献
默认情况下Endnote是将文献插入到文章最后面的。若需要将文献插入到各个章节后面,则在Endnote中设置,例如:点击“Edit”——“Output Styles”——“Edit BMC genomics”——“Sections”——选中“Create a bibliography for each section”——退出保存该格式为另外一个名字,然后使用这个保存的格式。

3.

human genome h38 infromation downloading

Writing date: 2015-11-17.

The latest Human Genome assembly version is : GRCh38 (GCA_000001405.15) . GRch38: Genome Reference Consortium Human Reference 38.

The GRch38 genome browses:
UCSC http://genome.ucsc.edu/cgi-bin/hgGateway
Ensembl http://www.ensembl.org/Homo_sapiens/Info/Index
Vega http://vega.sanger.ac.uk/Homo_sapiens/Info/Index
GENCODE http://www.gencodegenes.org/human_biodalliance.html

The downloading website of GRch38 information in Ensembl: http://www.ensembl.org/info/data/ftp/index.html
I recommend to download gh38 sequence functional annotations from Ensembl: ftp://ftp.ensembl.org/pub/release-82/genbank/homo_sapiens/.

mdkir sequence_annotation
cd sequence_annotation
lftp -e "mirror -c --parallel=5 /pub/release-82/genbank/homo_sapiens/" ftp://ftp.ensembl.org
cd ..

The downloading website of GRch38 information in GENCODE: http://www.gencodegenes.org/releases/23.html
I recommend to download gh38 fasta and gff3 files from GENCODE. These 2 files would be the main fasta and gff3 files for most users.

wget ftp://ftp.sanger.ac.uk/pub/gencode/Gencode_human/release_23/GRCh38.primary_assembly.genome.fa.gz
wget ftp://ftp.sanger.ac.uk/pub/gencode/Gencode_human/release_23/gencode.v23.basic.annotation.gff3.gz

SVG更改坐标系原点位置

在使用FigTree画树后。由于设置字体大小>14,于是导致export出来的图片中最上面一行字被截断了,从而使图片很丑。于是export出SVG格式文件。然后修改SVG坐标系原点位置,将图片完整显示出来。

在 <svg xmlns… 这行尾部添加 transform=”translate(0,20)” 解决。

纤维素,半纤维素和果胶的成份及其降解酶

1. 纤维素

Cellulose is a dominant structural polysaccharide in plants composed ofβ -D-glucose units with β-1,4-linkages.

Cellulose decomposition requires multiple enzymes. In general, cellulose is degraded to cellodextrin or cellobiose by the synergistic action of two cellulases: endoglucanase (EC 3.2.1.4) and cellobiohydrolase (EC 3.2.1.91) (Tomme et al., 1995; Bayer et al., 1998). Degradation of cellodextrin or cellobiose into monomeric glucose units requires another enzyme, β-glucosidase (EC 3.2.1.21), that hydrolyzes non-reducing 1,4-linked-β-glucose (Henrissat et al., 1989).

2. 半纤维素

Cellulose fibers are cross-linked by other polysaccharides called `hemicelluloses’ to increase the physical strength of the cell wall. Hemicelluloses include xylan (β-D-xylose units with β-1,4-linkages), glucomannan (β-D-mannose units andβ -D-glucose units with β-1,4-linkages), xyloglucan (β-D-glucose units with β-1,4-linkages, andβ -D-xylose and β-D-glucose units withβ -1,6-linkages), 1,3-1,4-β-glucan (β-D-glucose units with β-1,3- and β-1,4-linkages), and a relatively small amount of other polysaccharides composed of β-D-glucose,β -D-xylose, β-D-mannose and other sugar units with various linkages (McNeill et al., 1984).

3. 果胶

The scaffold of cellulose and hemicelluloses is filled with pectin (α-D-galacturonic acid units with mainly α-1,4-linkages), which functions as a cement-like substance in the cell wall.

reference:
Sakamoto, Kentaro, and Haruhiko Toyohara. “A comparative study of cellulase and hemicellulase activities of brackish water clam Corbicula japonica with those of other marine Veneroida bivalves.” Journal of Experimental Biology 212.17 (2009): 2812-2818.

通过WIG格式将转录组数据展示到Gbrowse2中

1. WIG格式介绍

WIG格式(Wiggle Track Format),可用于将转录组数据进行可视化展示。bigWig格式则是WIG格式的二进制方式,可以使用wigToBigWig将WIG格式转换成BigWig格式。
一个 WIG 格式实例文件:

track type=wiggle_0 name="sampleA1" description="RNA-Seq read counts of species A"
variableStep chrom=chr01 span=10
10001    13
10011    15
10021    12
fixedStep chrom=chr01 start=100031 step=10 span=10
17
15
20

如上例子,有2个注意点:

1. 第一行必须如理示例中格式。只有name和description这两个参数的值可以随意填写。
2. 有两种方法进行数据描述。分别是variableStep和fixedStep。前者数据内容用2行表示,后者数据部分仅用1行表示。
3. 这两种方法的几个参素意义为:
    chrom    设置序列名
    start    fixStep中Locus的起始位置
    step     fixStep中Locus的步进
    span     一个数据对应碱基数目

2. 将Bam文件转换成WIG文件并进行压缩

使用bam2wig命令将bam文件转换成wig文件。bam2wig命令可以来自于Augustus软件。

$ bam2wig sampleA1.tophat.bam > sampleA1.wig

该wig文件的span参数值为1。因此,当基因组越大,则wig文件越大。可以考虑设置span的值为10,能有效减小wig文件的大小。例如编写如些perl程序进行压缩wig文件。

#!/usr/bin/perl
use strict;

my $usage = <<USAGE;
Usage:
    perl $0 RNA-Seq.wig > RNA-Seq.cutdown.wig
USAGE
if (@ARGV==0){die $usage}

open IN, $ARGV[0] or die $!;

$_ = <>;
print;

my $locus = 1;
my $count = 0;
while () {
    if (m/^variableStep/) {
        $count = int(($count + 0.5) / 10);
        print "$locus\t$count\n" if $count > 0;
        s/$/ span=10/;
        print;
        $locus = 1;
    }
    else {
        if (m/(\d+)\s+(\d+)/) {
            my ($num1, $num2) = ($1, $2);
            if ($num1 >= $locus + 10) {
                $count = int(($count + 0.5) / 10);
                print "$locus $count\n" if $count > 0;
                $locus = $num1;
                $count = 0;
            }
            $count += $num2;
        }
    }
}

3. 将wig文件转换成wig binary文件和一个gff3文件

使用Gbrowse2所带命令 wiggle2gff3.pl 将wig文件转换成wig binary文件和一个gff3文件。每个基因组序列得到一个二进制格式的wig文件。同时生成一个gff3文件。该gff3文件指向所有的wig binary文件。

$ mkdir $PWD/gbrowse_track_of_RNA_seq
$ wiggle2gff3.pl --source=sampleA1 --method=RNA_Seq --path=$PWD/gbrowse_track_of_RNA_seq --trackname=track_A1 sampleA1.wig > sampleA1.gff3

4. 导入gff3文件到数据库,并配置Gbrowse配置文件

导入gff3文件

$ bp_seqfeature_load.pl -a DBI::mysql -d gbrowse2_species -u train -p 123456 sampleA1.gff3

配置文件:

[RNA_Seq_sampleA1_xyplot]
feature        = RNA_Seq:sampleA1
glyph          = wiggle_xyplot
graph_type     = boxes
height         = 50
scale          = right
description    = 1
category       = RNA-Seq:sampleA1
key            = Transcriptional Profile

[RNA_Seq_sampleA1_density]
feature        = RNA_Seq:sampleA1
glyph          = wiggle_density
height         = 30
bgcolor        = blue
description    = 1
category       = RNA-Seq:sampleA1
key            = Transcriptional Profile