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 | headLecture 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非常适合有生物学基础且想自学生信的这部分人入门。译者正是通过该课程入门的生信。
可点击 阅读原文 或在公众号的 菜单→课程系列 中找到本系列课程。
欢迎关注我们