1. SAMtools

简介SAMtools是用于处理SAM (Sequence Alignment/Map)格式文件的一系列工具,主要用来储存大容量的核酸测序结果。BAM是SAM文件的binary格式文件。SAMtools的主要作者是Heng Li,Heng Li在2012年因为对二代测序领域的贡献获得Benjamin Franklin Award

平台:Mac OS/Linux

安装:

  • 安装依赖的Ncurse库
1
# yum install ncurses*

快速运行

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
# sam格式文件转换为bam格式文件
# 新版本samtools不用使用-S
$ samtools view -b samFile > bamFile

# bam格式文件转换为sam格式文件
$ samtools view -h bamFile > samFile

# bam格式文件转换为sorted bam格式文件,用于长期储存和后续分析
# 后一个‘sortedBamFilePrefix’是指需要存储文件名前缀,比如想存储“human_1.bam”,则输入“human_1”
$ samtools sort bamFile sortedBamFilePrefix

# bam文件按照reads名称排序
$ samtools sort -n bamFile sortedBamFilePrefix

# 直接查看bam文件
$ samtools view bamFile | head -2

# 创建bam的index文件
$ samtools index bamFile

# 输出alignment数目,配合-f和-F过滤reads
$ samtools view -c bamFile

# 统计bam文件map数量,可以用于评估mapping的质量。需要输入indexed的bam文件
# 输出结果每一列分别为参考序列名称、参考序列长度、map上的reads数目、未map上的reads数目
$ samtools idxstats sortedIndexedBamFile

重要参数解释

  • view -h:打印bam文件头部,文件头部信息用于转换sam文件至bam文件。view -H:只打印bam文件头部。

补充

  • 一个sam文件式例子
A small sample of sam file
1
2
3
4
5
6
7
8
9
10
11
$ head samFile
HD      VN:1.0          SO:coordinate
@SQ     SN:chr1         LN:249250621
@SQ     SN:chr10        LN:135534747
@SQ     SN:chr11        LN:135006516

HISEQ2000-02:436:C2PG3ACXX:3:1316:4503:75198    99      chr1    11527   3       100M    =       11624   197     CTGGGTTTAAAAGTAAAAAATAAATATGTTTAATTTGTGAACTGATTACCATCAGAATTGTACTGTTCTGTATCCCACCAGCAATGTCTAGGAATGCCTG  ?@@FFFD?D>DDHFGGEGGIJDHJCEGCHE@FCHIJJJGG@GEGHGJGHHIIIJJJIJJECGG@FHIIJJJIIGEADHHEEDF;BB@CACCCEDDDDDAC    AS:i:0  XN:i:0  XM:i:0  XO:i:0  XG:i:0  NM:i:0  MD:Z:100     YT:Z:UU  NH:i:2  CC:Z:chr15      CP:i:102519544  HI:i:0

HISEQ2000-06:325:C2RC0ACXX:5:1101:16521:147635  355     chr1    11554   0       100M    =       11681   227     GTTTAATTTGTTAACTGATTACCATCAGAATTGTACTGTTCTGTATCCCACCAGCAATGTCTAGGAATGCCTGTTTCTCCACAAAGTGTTTACTTTTGGA  @@@DBDDDHFDDDHGGBA@AEFBBCGIC?BCGIAHBEFFHIHIDHC@DEGIIIHEG;8BDBFGGICFI:FG7@FHEHIICIEA;CA;77;7;;@;@A@6(    AS:i:-5 XN:i:0  XM:i:1  XO:i:0  XG:i:0  NM:i:1  MD:Z:11G88   YT:Z:UU  NH:i:5  CC:Z:chr12      CP:i:93956      HI:i:0

HISEQ2000-06:325:C2RC0ACXX:5:2304:4393:52082    163     chr1    11579   3       100M    =       11643   164     CAGAATTGTACTGTTCTGTATCCCACCAGCAATGTCTAGGAATACCTGTTTCTCCACAAAGTGTTTACTTTTGGATTTTTGCCAGTCTAACAGGTGAAGC  CCCFFFFFHHHHHIJJJJJJJJJJJJJJJJJJJJIJJJJJJJJJJJJJJJGIJJJIJJJJJFHFHHHIJJJJJIEHIJHHHHFFFFFFEEEEEDDDDDDD    AS:i:-6 XN:i:0  XM:i:1  XO:i:0  XG:i:0  NM:i:1  MD:Z:43G56   YT:Z:UU  NH:i:2  CC:Z:chr15      CP:i:102519492  HI:i:0

2. Bowtie2

Bowtie使用介绍,详见二代测序中的短序列比对

3. Trinity

TrinityBroad Institute开发的根据RNA-seq数据从头(de novo)组装转录组的工具。

4. miRDeep2

miRDeep2 是用于高通量测序中检测miRNA 的工具。

5. MISO

简介MISO(Mixture-of-Isoforms)是用来计算和探测RNA-seq数据中不同样本的基因选择剪切。

平台:Python跨平台使用。

快速运行

An Example of Running MISO
1

6. 质量检测

6.1 RNA-SeQC

RNA-SeQC是用于检测RNA-seq测序质量。

6.2 RSeQC

RSeQC用于对RNA-seq数据质量控制。

7. 小工具集锦

7.1 操作fastq文件

在Linux机器上配以一系列命令,操作fastq或者fastq.gz文件,比如mySeq_1_1.fastq测序文件。

Manipulating fastq file with bash
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
# Fastq文件前几行
$ head mySeq_1_1.fastq
@HISEQ2000-06:325:C2RC0ACXX:5:1101:1217:2215 1:N:0:TGCCGGCT
GGCACAGTCCATGCTTTTAACCAGATTTGAACAGAAGAATGGCCACTTGNNNNNNGTAGNNNNNNANNNNGNNNTNNNTNNCATGTGTCACATAACTACC
+
@C@FFFFDFFHHHGFFIIBEGHJIFHIIGCGIJJIJIJJIJIJJIIIGI######-.<B######-####,###,###(##,,5<?CDDDDCCDDDCDDD
@HISEQ2000-06:325:C2RC0ACXX:5:1101:1621:2040 1:N:0:TGCCGGCT
NGGTGACCTTCTCCCGCCAGAAGCCAAAAACTGCAGCCTACTTTTCTGAAGTGGTTATCTTGGGACTGAGGTATGGGCTATCTTGGGCTGTTTCCTATTT
+
#1==D;ADBHDFDHFHGGIIB;FFHIIGGGIGADEGEIIAHEEHIGHH==C8BGC;=FCDHICHCAAC>CE..;A>AC?9@CCCCCBBCB?CCCCCCCCC
@HISEQ2000-06:325:C2RC0ACXX:5:1101:1663:2146 1:N:0:TGCCGGCT
AGAGCCAAATATTTCAACAAAACTGCAGTTTAATTTCAGAAAATGTTAAAATATATATTTATACATCAATTTCTGACATACACTTAATGTGTTAGTATAC

# 统计fastq文件中的reads数量
$ awk 'END{print NR/4}' mySeq_1_1.fastq

# 查看fastq.gz文件
$ zcat mySeq_1_1.fastq.gz

# 统计fastq.gz文件reads数量
# 以下两种方法等价
$ zcat mySeq_1_1.fastq.gz | awk 'END{print NR/4}'
$ zcat mySeq_1_1.fastq.gz | echo $((`wc -l`/4))
# 同时,由于一些测序文件会对read做特殊标记,比如“@HISEQ2000”。所以,可以使用正则匹配。
# 注意:这种方法有潜在错误,慎重使用。
$ zgrep -c '@HISEQ200' mySeq_1_1.fastq

# 统计fastq文件reads长度分布
$ zcat mySeq_1_1.fastq.gz | awk '{if(NR%4==2) print length($1)}' | sort -n | uniq -c

7.2 操作fasta文件

Manipulating fasta file with bash
1
2
3
4
5
6
# 序列数目
$ grep -c '>@' mySeq.fa

# 提取特定序列
samtools faidx mySeq.fa
samtools faidx mySeq.fa chr1 chr2 chr3

7.3 Picard

Picard是一个Java平台的工具包,包括一系列处理高通量测序的命令行工具。

7.4 FASTX-Toolkit

FASTX-Toolkit是一系列用于处理大量短序列FASTA和FASTQ文件的工具。

8. 各种数据类型

为了下游分析,高通量测序结果往往使用多种数据格式储存,详细参考File Formats

8.1 GTF和GFF文件互转

使用Cufflinksgffread直接进行转换。

Convert the File Format Between GTF and GFF
1
2
3
4
5
# GTF转换为GFF/GFF3
$ gffreads myGtfFile.gtf -o myGffFile.gff3

# GFF/GFF3转换为GTF
$ gffread myGffFile.gff3 -T -o myGtfFile.gtf

8.2 GTF文件提取序列信息

使用Cufflinksgffread直接进行转换。

Extract sequences from GTF files
1
2
# 从GTF文件中提取序列信息
$ gffread -w myFastaFile.fa -g myGenome.fa myGtfFile.gtf

参考资料

更新记录

2018年4月3日

Comments