生物信息 BIOINFORMATICS Comments turned off

在Ubuntu环境中设置Bowtie

By on 8月 20th, 2012

什么是Bowtie？What is Bowtie？

Bowtie is an ultrafast, memory-efficient short read aligner geared toward quickly aligning large sets of short DNA sequences (reads) to large genomes. It aligns 35-base-pair reads to the human genome at a rate of 25 million reads per hour on a typical workstation. Bowtie indexes the genome with a Burrows-Wheeler index to keep its memory footprint small: for the human genome, the index is typically about 2.2 GB (for unpaired alignment) or 2.9 GB (for paired-end or colorspace alignment). Multiple processors can be used simultaneously to achieve greater alignment speed. Bowtie can also output alignments in the standard SAM format, allowing Bowtie to interoperate with other tools supporting SAM, including the SAMtools consensus, SNP, and indel callers. Bowtie runs on the command line under Windows, Mac OS X, Linux, and Solaris.

Bowtie also forms the basis for other tools, including TopHat: a fast splice junction mapper for RNA-seq reads, Cufflinks: a tool for transcriptome assembly and isoform quantitiation from RNA-seq reads, Crossbow: a cloud-computing software tool for large-scale resequencing data,and Myrna: a cloud computing tool for calculating differential gene expression in large RNA-seq datasets.

If you use Bowtie for your published research, please cite the Bowtie paper.
简而言之，Bowtie是分析比对高通量（下一代）测序技术数据的一个工具软件。
Bowtie的安装
[bash]sudo apt-get install bowtie[/bash]
设置Bowtie的Index所在位置
[bash]gedit ~/.bashrc[/bash]
在其中添加一行（请根据实际情况修改文件夹名称）：
[bash]export BOWTIE_INDEXES=”~/bowtie-index”;[/bash]
保存、关闭gedit。
重新读取设置：[bash]source ~/.bashrc[/bash]
这时运行：[bash]export|grep BOWTIE[/bash]
会得到BOWTIE_INDEXES的变量设置情况。以后既可以将Bowtie运行所需要的Index文件下载下来（或者自己生成）放置到该文件夹中。
如此这般，既可以使用Bowtie了。如果使用Bowtie时发生内存不够的错误，多加点内存吧，别无他法。
生物信息 BIOINFORMATICS, 生物学 BIOLOGY, 转载 SHARE Comments turned off

生物信息学中的 SAM 格式介绍

By on 8月 16th, 2012

Following is a brief description of the SAM format as output by bowtie when the -S/–sam option is specified. For more details, see the SAM format specification.

When -S/–sam is specified, bowtie prints a SAM header with @HD, @SQ and @PG lines. When one or more –sam-RG arguments are specified, bowtie will also print an @RG line that includes all user-specified –sam-RG tokens separated by tabs.

Each subsequnt line corresponds to a read or an alignment. Each line is a collection of at least 12 fields separated by tabs; from left to right, the fields are:

  1. Name of read that aligned
  2. Sum of all applicable flags. Flags relevant to Bowtie are:
    1 The read is one of a pair
    2 The alignment is one end of a proper paired-end alignment
    4 The read has no reported alignments
    8 The read is one of a pair and has no reported alignments
    16 The alignment is to the reverse reference strand
    32 The other mate in the paired-end alignment is aligned to the reverse reference strand
    64 The read is the first (#1) mate in a pair
    128 The read is the second (#2) mate in a pair

    Thus, an unpaired read that aligns to the reverse reference strand will have flag 16. A paired-end read that aligns and is the first mate in the pair will have flag 83 (= 64 + 16 + 2 + 1).

  3. Name of reference sequence where alignment occurs, or ordinal ID if no name was provided
  4. 1-based offset into the forward reference strand where leftmost character of the alignment occurs
  5. Mapping quality
  6. CIGAR string representation of alignment
  7. Name of reference sequence where mate’s alignment occurs. Set to = if the mate’s reference sequence is the same as this alignment’s, or * if there is no mate.
  8. 1-based offset into the forward reference strand where leftmost character of the mate’s alignment occurs. Offset is 0 if there is no mate.
  9. Inferred insert size. Size is negative if the mate’s alignment occurs upstream of this alignment. Size is 0 if there is no mate.
  10. Read sequence (reverse-complemented if aligned to the reverse strand)
  11. ASCII-encoded read qualities (reverse-complemented if the read aligned to the reverse strand). The encoded quality values are on the Phred quality scale and the encoding is ASCII-offset by 33 (ASCII char !), similarly to a FASTQ file.
  12. Optional fields. Fields are tab-separated. For descriptions of all possible optional fields, see the SAM format specification. bowtieoutputs some of these optional fields for each alignment, depending on the type of the alignment:
  • NM:i:<N>    Aligned read has an edit distance of <N>.
  • CM:i:<N>    Aligned read has an edit distance of <N> in colorspace. This field is present in addition to the NM field in -C/–color mode, but is omitted otherwise.
  • MD:Z:<S>    For aligned reads, <S> is a string representation of the mismatched reference bases in the alignment. See SAM format specification for details. For colorspace alignments, <S> describes the decoded nucleotide alignment, not the colorspace alignment.
  • XA:i:<N>    Aligned read belongs to stratum <N>. See Strata for definition.
  • XM:i:<N>    For a read with no reported alignments, <N> is 0 if the read had no alignments. If -m was specified and the read’s alignments were supressed because the -m ceiling was exceeded, <N> equals the -m ceiling + 1, to indicate that there were at least that many valid alignments (but all were suppressed). In -M mode, if the alignment was randomly selected because the -M ceiling was exceeded, <N> equals the -M ceiling + 1, to indicate that there were at least that many valid alignments (of which one was reported at random).

 

SAM文件定义

SAM是一种序列比对格式标准， 由sanger制定，是以TAB为分割符的文本格式。主要应用于测序序列mapping到基因组上的结果表示，当然也可以表示任意的多重比对结果。

不同的软件，不同的时期，不同的研究方向，都会创建一种或者多种格式标准，当然根据当时的需要，创建符合需求的标准，也是最容易的事情，而反过来想 要真正的理解标准，也必须理解为什么要创建这样的标准，解决什么样的需要。我前面的有篇文章已经对于现有的多重比对的格式进行总结，但其更多的站在比较基 因组学的角度。当我们去了解sam标准格式是什么的时候，就要思考既然以及有了这么多得标准，为什么还要定义SAM标准，当然拿所有的格式进行比较也并非 易事，但是简单的对比，就可以了解其中一二，比如aln格式，是比对视图化的展示，存储的信息不够结构化，无法方便的作为另外程序的输入；表示信息的有限 性，如果100个多重比对序列放到一个文件中，查看维护就会非常困难；还有些格式标准挺强大，但是太繁琐，同时不够灵活。那么反过来就是SAM格式的优 点，那么SAM如何做到这一点的呢？

SAM要处理好的问题：

  • 非常多序列（read)，mapping到多个参考基因组（reference）上；
  • 同一条序列，分多段（segment）比对到参考基因组上；
  • 无限量的，结构化信息表示，包括错配、删除、插入等比对信息；

SAM分为两部分，注释信息（header section）和比对结果部分（alignment section），注释信息可有可无，都是以@开头，用不同的tag表示不同的信息，主要有@HD，说明符合标准的版本、对比序列的排列顺序；@SQ，参 考序列说明；@RG，比对上的序列（read）说明；@PG，使用的程序说明；@CO，任意的说明信息。

SAM格式示例

比对结果部分（alignment section），每一行表示一个片段（segment）的比对信息，包括11个必须的字段（mandatory fields）和一个可选的字段，字段之间用tag分割。必须的字段有11个，顺序固定，不可用时，根据字段定义，可以为’0‘或者’*‘，这是11个字 段包括：

  1. QNAME，比对片段的（template）的编号；
  2. FLAG，位标识，template mapping情况的数字表示，每一个数字代表一种比对情况，这里的值是符合情况的数字相加总和；
  3. RNAME，参考序列的编号，如果注释中对SQ-SN进行了定义，这里必须和其保持一致，另外对于没有mapping上的序列，这里是’*‘；
  4. POS，比对上的位置，注意是从1开始计数，没有比对上，此处为0；
  5. MAPQ，mappint的质量；
  6. CIGAR，简要比对信息表达式（Compact Idiosyncratic Gapped Alignment Report），其以参考序列为基础，使用数字加字母表示比对结果，比如3S6M1P1I4M，前三个碱基被剪切去除了，然后6个比对上了，然后打开了一 个缺口，有一个碱基插入，最后是4个比对上了，是按照顺序的；
  7. RNEXT，下一个片段比对上的参考序列的编号，没有另外的片段，这里是’*‘，同一个片段，用’=‘；
  8. PNEXT，下一个片段比对上的位置，如果不可用，此处为0；
  9. TLEN，Template的长度，最左边得为正，最右边的为负，中间的不用定义正负，不分区段（single-segment)的比对上，或者不可用时，此处为0；
  10. SEQ，序列片段的序列信息，如果不存储此类信息，此处为’*‘，注意CIGAR中M/I/S/=/X对应数字的和要等于序列长度；
  11. QUAL，序列的质量信息，格式同FASTQ一样。

可选字段（optional fields)，格式如：TAG:TYPE:VALUE，其中TAG有两个大写字母组成，每个TAG代表一类信息，每一行一个TAG只能出现一次，TYPE表示TAG对应值的类型，可以是字符串、整数、字节、数组等。

要注意的几个概念，以及与之对应的模型：

  • reference
  • read
  • segment
  • template（参考序列和比对上的序列共同组成的序列为template）
  • alignment
  • seq

 

http://bowtie-bio.sourceforge.net/manual.shtml#sam-bowtie-output

http://boyun.sh.cn/bio/?p=1890

信息技术 INFORMATICS, 原创教程 TUTORIAL, 矩阵实验室 MATLAB Comments turned off

MatLab中国地图进阶：以不同的颜色表示各省数据的差异

By on 6月 5th, 2012

使用MatLab绘制中国地图十分方便，我在前文中已经有相应阐述。实际应用中，我们常常需要将不同的区域以不同的颜色展示出来，这就是本文要阐述的内容。

本文使用MatLab将各省人口数据展示出来，人口越多颜色越黄，人口越少颜色越绿。最终结果如下图所示（点击看大图）。

绘制过程如下：

%-- 分省中国地图在MATLAB中的画法 --%
% Author: gaospecial@gmail.com
% Date: 2012-6-5 14:04:01
% 地图数据来源：国家基础地理信息系统网站
% 人口数据来源：国家统计年鉴2011版

载入地图数据

载入各省的多边形数据

sheng=shaperead('maps/bou2_4p.shp', 'UseGeoCoords', true); % 省
% 检查数据发现中午名称为乱码，可能是由于字符编码引起的，需要纠正文件导入时出现的乱码
% 中文名称记录在相应的dbf文件中，可以使用Excel打开，并且不会出现乱码，截取其中的中文保存为文本文件
% 然后在导入工作区，并将其保存为*.mat格式文件，随用随取。
% 依次替换sheng数据中的乱码【算是变通之法，有谁知道更好的方法一定告诉我哦】
load chinese_name.mat% 省，省会，主要城市的正确中文
for i=1:length(sheng)
    sheng(i).NAME=sheng_chinese_name{i}; % 纠正中文显示错误
end
unique(sheng_chinese_name) % 含有34个省（直辖市）的数据
length(sheng) % 共分为925个区块
ans = 

    '上海市'
    '云南省'
    '内蒙古自治区'
    '北京市'
    '台湾省'
    '吉林省'
    '四川省'
    '天津市'
    '宁夏回族自治区'
    '安徽省'
    '山东省'
    '山西省'
    '广东省'
    '广西壮族自治区'
    '新疆维吾尔自治区'
    '江苏省'
    '江西省'
    '河北省'
    '河南省'
    '浙江省'
    '海南省'
    '湖北省'
    '湖南省'
    '甘肃省'
    '福建省'
    '空'
    '西藏自治区'
    '贵州省'
    '辽宁省'
    '重庆市'
    '陕西省'
    '青海省'
    '香港特别行政区'
    '黑龙江省'

ans =

   925

导入2011年各省人口数据

使用importdata向导导入2011年全国31个省的人口数据：data为人口数目 textdata为省名称

d=importdata('renkou.txt');
data=d.data; % 人口数目
textdata=d.textdata; % 相对应的省的名称

定义地图参数

针对不同省份，分别设置不同的颜色（FaceColor） 定义颜色

k=128;
mycolormap=summer(k);
% 生成不同区域按大小的颜色，按照人口数目多少分别指定不同的颜色
% 人口越多，颜色越突出
geoname={sheng.NAME}';
max_data = max(data);
n=length(data);
mysymbolspec=cell(1,n);% 预定义变量可以加快处理速度
for i=1:n
    count=data(i);
    mycoloridx=floor( k * count / max_data );
    mycoloridx(mycoloridx<1)=1;
    myprovince=textdata{i};
    geoidx=strmatch(myprovince, geoname);
    if numel(geoidx) > 0
        province_name=geoname( geoidx(1) );
        mysymbolspec{i} = {'NAME', char(province_name), 'FaceColor', mycolormap( mycoloridx, :) };
    end
end

显示地图

figure
ax=worldmap('china'); % 使用worldmap的坐标轴作图
setm(ax,'grid','off') % 关闭grid
setm(ax,'frame','off') % 关闭边框
setm(ax,'parallellabel','off') % 关闭坐标轴标记
setm(ax,'meridianlabel','off') % 关闭坐标轴标记

% 最关键的两个语句
symbols=makesymbolspec('Polygon',{'default','FaceColor',[0.9 0.9 0.8],...
    'LineStyle,'--','LineWidth',0.2,...
    'EdgeColor',[0.8 0.9 0.9]},...
    mysymbolspec{:}...
    );
geoshow(sheng,'SymbolSpec',symbols); % 此处用mapshow投影会不正确

图的标注

% 在图像右侧显示bar
colormap(summer(k))
hcb=colorbar('EastOutside');
step=round(max_data/11);
set(hcb,'YTick',(0:.1:1))
set(hcb,'YTickLabel',num2cell(0:step:max_data))

% 给图像加标题
title('中国2011年各省人口数（万人）')

下载本文附件，包括M文件和输出结果（附件已更新2012-6-14）。
china_map

信息技术 INFORMATICS, 原创教程 TUTORIAL, 网站建设 WEBSITE Comments turned off

the_content()和the_excerpt()

By on 6月 5th, 2012

这是wordpress中用以提取post内容的常用函数。它们用法基本相同，很多时候可以相互取代，分别用于显示post全部内容和post摘要。

然而，事实上在应用中它们最大的区别在于，前者输出格式化的内容，而后者则输出文字（格式错乱，没有超链接，换行，样式定义等）。所以，大多数时候，我只用前一个函数。

有些主题作者喜欢用the_excerpt()（特别是在分类列表中），我常常将其改为the_content()。

the_content   Display the post content.
void the_content ([string $more_link_text = null], [bool $stripteaser = false])
string $more_link_text: Optional. Content for when there is more text.
bool $stripteaser: Optional. Strip teaser content before the more text. Default is false.

the_excerpt  Display the post excerpt.
uses: apply_filters() – Calls ‘the_excerpt’ hook on post excerpt.
void the_excerpt ()

 

信息技术 INFORMATICS, 原创教程 TUTORIAL, 生物信息 BIOINFORMATICS, 生物软件 BIO-SOFTWARE, 自然科学 NATURE & SCIENCE, 转载 SHARE Comments turned off

删除Zotero自动生成的Tags或关键词

By on 3月 14th, 2012

添加项目时，Zotero默认将关键词以Tag的形式加入，这样可以快速的查看与某一关键词相关的条目。可是自动加入的Tag或关键词有太多太多，时间一长，Zotero数据库中会产生成百上千个关键词或Tag，导致其失去了便于检索的应用价值。

所以，建议在取消该选项，而改为人工添加更加简单易记的关键词或Tag。

那么，对于已经存在的Tag如何一下子全部删除呢？按照下面的步骤操作即可。

  • 打开Zotero配置页面，选择高级（Advanced），点击“打开数据名录”（Show data directory）；
  • 关闭Firefox，关闭之前将本文用IE打开，或复制到一个文档中，以便随时参考；
  • 复制“数据目录”中的“zotero.sqlite”为“zotero.sqlite.bak”，即备份Zotero数据库；
  • 下面对Zotero数据库进行查询，可以使用任何熟悉的数据库工具，这里以“SharpPlus SQlite Developer”，操作如下：
  1. 注册数据库，根据向导选择zotero.sqlite文件。
  2. 双击导入的数据库。
  3. 选择“工具”-“SQL编辑器”，在其中输入以下代码：
    [cc lang=”sql”]delete from itemTags where itemTags.tagID in (select
    tags.TagID from tags where type = 1);
    delete from tags where type = 1;[/cc]
  4. 选择“运行”，即窗口中第一个图标。
  5. 这时候，窗口下半部分应当有信息出现“Query OK…”等等。
  • 关闭该工具，重新打开Firefox。
  • 完成。

Refer：http://forums.zotero.org/discussion/4051/remove-all-tags/