Biostar（大结局）：课程29、30

翻译：生物女博士

注：本系列课程翻译已获得原作者授权。

# 为作者的注释

#* 为译者的注释

Lecture 29 - 比较特征计数软件

# 下载和安装 eXpress

cd ~/src

#* Mac

curl -OL http://bio.math.berkeley.edu/eXpress/downloads/express-1.5.1/express-1.5.1-macosx_x86_64.tgz
tar zxvf express-1.5.1-macosx_x86_64.tgz
ln -fs ~/src/express-1.5.1-macosx_x86_64/express ~/bin

# Linux

# curl -OL http://bio.math.berkeley.edu/eXpress/downloads/express-1.5.1/express-1.5.1-linux_x86_64.tgz
#* tar zxvf express-1.5.1-linux_x86_64.tgz
#* ln -fs ~/src/express-1.5.1-linux_x86_64/express ~/bin

# 来看一下eXpress程序是否能运行

express

# 安装subread包

cd ~/src
curl -OL http://downloads.sourceforge.net/project/subread/subread-1.4.6/subread-1.4.6-source.tar.gz
tar zxvf subread-1.4.6-source.tar.gz
cd subread-1.4.6-source/src
make -f Makefile.MacOS #* 如果是linux系统，则改为 make -f Makefile.Linux

#* 根据自己的系统来修改

ln -fs ~/src/subread-1.4.6-source/bin/featureCounts ~/bin ln -fs ~/src/subread-1.4.6-source/bin/subread-align ~/bin ln -fs ~/src/subread-1.4.6-source/bin/subread-buildindex ~/bin

# 看一下featureCounts程序能否正常运行

featureCounts
subread-align

# 数据准备。

# 假如你还没下边的数据，现在先来获取数据。

curl -OL http://www.personal.psu.edu/iua1/courses/files/2014/rnaseq-data.tar.gz
tar zxvf rnaseq-data.tar.gz
rm -f rnaseq-data.tar.gz

# 如果你没有GTF文件，先下载一下。

curl -L ftp://ftp.ensembl.org/pub/release-77/gtf/homo_sapiens/Homo_sapiens.GRCh38.77.gtf.gz | gunzip -c > ~/refs/hg38/hg38.gtf.gz

# 只提取GTF文件里的4号染色体信息。

# 因为我们知道我们的比对数据只是4号染色体的。

cat ~/refs/hg38/hg38.gtf.gz | awk ' $1=="4" { print $0 } ' > chr4.gtf

# 使用TopHat定量

#* 这里用到的人类基因组4号染色体序列~/refs/hg38/chr4是在上节课下载、解压和索引过的，详情请参照课程27和28。

tophat -G chr4.gtf -o data/ctl1-tophat ~/refs/hg38/chr4 data/ctl1.fastq.gz
tophat -G chr4.gtf -o data/ctl2-tophat ~/refs/hg38/chr4 data/ctl2.fastq.gz

# 对比这两个对照样品……看看会发生什么。是否应该有表达差异？

cuffdiff -o diff chr4.gtf data/ctl1-tophat/accepted_hits.bam data/ctl2-tophat/accepted_hits.bam

# 使用eXpress定量

# eXpress aligns against a transcript.

# 我们可以从数据源下载转录本数据，或者从已有的注释中生成。

# CuffLinks的包中含有一个GFF FASTA格式软件

# 这需要用到参考FASTA序列和chr4.gtf文件

ln -sf ~/src/cufflinks-2.1.1.OSX_x86_64/gffread ~/bin

# 运行查看帮助文档

gffread -h

# 根据外显子坐标提取序列。

gffread chr4.gtf -t exon -g ~/refs/hg38/chr4.fa -w ~/refs/hg38/chr4-transcripts.fa

# 我们也可以从Ensembl上直接获取转录本

# ftp://ftp.ensembl.org/pub/release-77/fasta/homo_sapiens/cdna/Homo_sapiens.GRCh38.cdna.all.fa.gz
# gunzip Homo_sapiens.GRCh38.cdna.all.fa.gz

# eXpress 需要使用一个比对文件来进行定量。

# 而这个比对文件需要以read名字排序。

# 通过bwa创建一个比对文件。你也可以使用别的比对软件。

mkdir -p bwa
bwa index ~/refs/hg38/chr4-transcripts.fa

# 注意这个BAM文件不是以位置排序的。

# 这将无法在IGV进行可视化！

bwa mem ~/refs/hg38/chr4-transcripts.fa data/ctl1.fastq.gz | samtools view -b - > bwa/ctl1.bam

# 用 express 进行定量

express ~/refs/hg38/chr4-transcripts.fa bwa/ctl1.bam -o bwa

# 用SubRead进行定量

# 这需要我们跟基因组比对！

# 用tophat创建一个比对。

tophat -G chr4.gtf -o data/ctl1-tophat ~/refs/hg38/chr4 data/ctl1.fastq.gz

# 用featureCounts统计特征（features）

mkdir -p sub
featureCounts -a chr4.gtf -g transcript_id data/ctl1-tophat/accepted_hits.bam  -o sub/ctl1.tophat.txt

# 我们来比较一下结果。

# featureCounts的统计值是基于TopHat比对

cat sub/ctl1.tophat.txt | sort -rn -k 7,7 | cut -f 1,6,7 | head | cat -n

# eXpress的统计值是基于bwa比对

cat bwa/results.xprs | sort -rn -k5,5 | cut -f 2,3,5 | head | cat -n

# 要从 featureCounts的输出结果中获取FPKM，你需要除以长度，并以文库大小进行标准化：60683

cat sub/ctl1.tophat.txt | grep "ENST" |  awk '{ print $1, $7/$6/60683 * 10^9 } ' | sort -rn -k2 | head

Lecture 30 - 基因本体论（Gene Ontology)

# 基因本体。了解数据的含义。

# 基于TopHat结果进行分析

# 有哪些列？

cat gene_exp.diff | head -1

# 我们有多种选择，根据不同列进行排序。

# 最大/最小的倍数改变（ fold change），样本1 或者样本2 的覆盖度等等

# 通常我们以p-value大小来处理基因

# 有时候你需要根据一些特性来分类基因：上调的、下调的以及对每个基因集分别进行富集分析。

# 倍数变化（foldchange）为正的基因

cat gene_exp.diff |  grep yes | awk '$10 > 0 { print $0 }' | cut -f 3,13 > positive.txt

# 倍数变化（foldchange）为负的基因

cat gene_exp.diff |  grep yes | awk '$10 < 0 { print $0 }' | cut -f 3,13 > negative.txt

# 也有一些工具，需要用整个数据集来计算所有基因数值的背景分布。

# 此时，我们取第2列 (ensemble ID).

cat gene_exp.diff | cut -f 2,13 > all-ensemble.txt

# 用这些列表，粘贴到GO ontology搜索的程序中。

# 安装 ErmineJ

cd ~/src
curl -OL http://www.chibi.ubc.ca/ermineJ/distributions/ermineJ-3.0.2-generic-bundle.zip
unzip ermineJ-3.0.2-generic-bundle.zip
ln -sf ~/src/ermineJ-3.0.2/bin/ermineJ.sh ~/bin

# 似乎，我们需要使它变得可执行。

chmod +x ~/bin/ermineJ.sh

# 用图形用户界面运行。

ermineJ.sh --gui

# ermine

Biostar非常适合有生物学基础且想自学生信的这部分人入门。译者正是通过该课程入门的生信。

可点击 阅读原文 或在公众号的 菜单→课程系列 中找到本系列课程。