搜索
查看: 2013|回复: 0

[Other] Biostar(大结局):课程29、30差异表达分析以及用ErmineJ做GO

[复制链接]

33

主题

46

帖子

230

积分

管理员

Rank: 9Rank: 9Rank: 9

积分
230
发表于 2017-9-24 14:47:58 | 显示全部楼层 |阅读模式
注:本文是生信媛微信公众号原创文章
作者:István Albert
原文链接:https://mp.weixin.qq.com/s/QMIfwPNPo6uCebvf1EOyIQ


翻译:生物女博士
注:本系列课程翻译已获得原作者授权。

# 为作者的注释
#* 为译者的注释

Lecture 29 - 比较特征计数软件

# 下载和安装 eXpress

cd ~/src

#* Mac
curl -OL http://bio.math.berkeley.edu/eXp ... acosx_x86_64.tgztar 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/eXp ... 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 ... .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 ... aseq-data.tar.gz
tar zxvf rnaseq-data.tar.gz
rm -f rnaseq-data.tar.gz

# 如果你没有GTF文件,先下载一下。
curl -L ftp://ftp.ensembl.org/pub/releas ... ns.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/releas ... h38.cdna.all.fa.gz
# gunzip Homo_sapiens.GRCh38.cdna.all.fa.gz

# eXpress 需要使用一个比对文件来进行定量。
# 而这个比对文件需要以read名字排序。
# 通过bwa创建一个比对文件。你也可以使用别的比对软件。
mkdir -p bwabwa 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/ ... ric-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非常适合有生物学基础且想自学生信的这部分人入门。译者正是通过该课程入门的生信。




欢迎到微信公众号订阅我们
生信媛
bio_sxy


本帖子中包含更多资源

您需要 登录 才可以下载或查看,没有帐号?立即注册

x



上一篇:Biostar:课程27、28:RNAseq分析
下一篇:如果一个新的ubuntu一定要安装一些库
回复

使用道具 举报

您需要登录后才可以回帖 登录 | 立即注册

本版积分规则

QQ|手机版|小黑屋|生信技能树 ( 粤ICP备15016384号  

GMT+8, 2019-7-21 04:55 , Processed in 0.033640 second(s), 27 queries .

Powered by Discuz! X3.2

© 2001-2013 Comsenz Inc.