搜索
查看: 2708|回复: 0

[Other] Biostar:课程27、28:RNAseq分析

[复制链接]

33

主题

46

帖子

230

积分

管理员

Rank: 9Rank: 9Rank: 9

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


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

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

Lecture 27 - RNA-Seq入门
#* 由于课程时间比较上比较旧,所以使用的还是TopHat
#* 目前已经用Hisat2取代了TopHat流程了。具体使用请参考我们公众号文章《新司机带你学RNA-Seq数据分析》by 徐春晖
#* 另外,使用STAR进行RNA-seq分析见《比对软件STAR的使用—高通量测序数据处理学习记录(一)》by 徐春晖
#* 更多RNA-seq分析相关文章可以使用我们公众号全文搜索程序(https://xuzhougeng.shinyapps.io/biosxy/)(洲更同学荣誉出品)进行搜索。

# 安装TopHat
cd ~/src

# Mac OSX

# Linux


# 解压和安装
tar zxvf tophat-2.0.13.OSX_x86_64.tar.gz
ln -s ~/src/tophat-2.0.13.OSX_x86_64/tophat ~/bin

# 跟着Stephen Turner的RNA-Seq训练营来学习。
#*  网址已经失效了。找到一个RNAseq分析相关的:http://bioconnector.org/workshops/r-rnaseq-airway.html

# 从Ensembl下载基因组数据。
# 为人类基因组创建一个目录
mkdir -p ~/refs/hg38

# 获取基因组序列

# 获取注释信息

# 解压hg38目录下的文件
gunzip ~/refs/hg38/*.gz

# 我们可以建一个简化版的只含4号染色体数据的GTF文件
cat ~/refs/hg38/hg38.gtf | awk ' $1 =='4' { print $0 }' > chr4.gtf

# 你还可以只提取注释文件的特定部分,比如外显子。
cat chr4.gtf | awk '$3=="exon" { print $0 } ' > exons.gtf

# 建立bowtie2版的参考序列。
bowtie2-build ~/refs/hg38/chr4.fa  ~/refs/hg38/chr4

# 建立bwa版的参考序列。
bwa index ~/refs/hg38/chr4.fa

# 获取和解压数据。

# 首先,我们来比对这个数据
# 给从sam到bam文件中的转换和排序代码创建一个缩写。
alias tobam='samtools view -b - | samtools sort -o - tmp '
time bwa mem ~/refs/hg38/chr4.fa  data/ctl1.fastq.gz | tobam > ctl1.bwa.bam
samtools index ctl1.bwa.bam


# 查看覆盖度最高的区域
bedtools coverage -abam ctl1.bwa.bam -b chr4.gtf | more

# 把它放到一个文件中。
bedtools coverage -abam ctl1.bwa.bam -b chr4.gtf > coverages.txt

# 以第10列排序(这一列的数据是map上的reads的数目)
# Tab分隔的文件需要像下边这样才能正常运行
# 否则sort会在任何的空格符处进行分开。
cat coverages.txt | sort -t$'\t' -rn -k10,10 | head

# 把数据导入IGF
# 找到基因RPS3A
# 我们使用相同的数据,用可以识别序列剪切位置的比对软件:TopHat
tophat ~/refs/hg38/chr4 data/ctl1.fastq.gz

# 索引一下结果文件
samtools index tophat_out/accepted_hits.bam

#到IGV查看结果。


Lecture 28 - 运行Tuxedo suite
#* Tuxedo suite主要包含了Bowtie、TopHat和Cufflinks
#* Bowtie建立索引:bowtie2-build ~/refs/hg38/chr4.fa  ~/refs/hg38/chr4
#* TopHat将RNA-seq的reads比对到参考序列上:tophat ~/refs/hg38/chr4
data/ctl1.fastq.gz

#* Cufflinks计算差异表达水平
# 获取和安装cuffdiff以及cufflinks
cd ~/src

# Mac版本的下载和安装
curl -OL http://cole-trapnell-lab.github. ... SX_x86_64.tar.gz
tar zxvf cufflinks-2.1.1.OSX_x86_64.tar.gz
ln -fs ~/src/cufflinks-2.1.1.OSX_x86_64/cufflinks ~/bin
ln -fs ~/src/cufflinks-2.1.1.OSX_x86_64/cuffdiff ~/bin


# Linux下则下载另一个不同版本。注意目录名字的改变。
# 自行根据上边的命令进行修改。

# 来看一下是不是可以正常使用了。
cuffdiff

# 我们需要用tophat来比对,然后用cuffdiff来进行定量。
# 我们假设基因组已经如lecture27那样,下载并且索引了
#* 这一步获取了一些rnaseq的数据
curl -OL http://www.personal.psu.edu/iua1 ... aseq-data.tar.gz
tar zxvf rnaseq-data.tar.gz
rm rnaseq-data.tar.gz


# 如果没有,那就通过下边链接获取。
# gunzip ~/refs/hg38/*.gz

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

# 对每个文件运行tophat
tophat -G chr4.gtf -o data/ctl1-tophat ~/refs/hg38/chr4 data/ctl1.fastq.gz

# 复制结果文件到bam目录,并换一个名字。
cp data/ctl1-tophat/accepted_hits.bam bam/ctl1.bam
samtools index bam/ctl1.bam


# 现在来重复之前的每个步骤
# 你应该写一个脚本,或者将所有命令列入一个更自动的脚本(见下文)
# bash run-tophat.sh

# 现在我们用cuffdiff来看一下差异表达。
cuffdiff

# 运行cuffdiff,并把运行结果放到bam/cuffdiff_out文件夹下
cd bam
cuffdiff -o cuffdiff_out chr4.gtf ctl1.bam,ctl2.bam,ctl3.bam uvb1.bam,uvb2.bam,uvb3.bam

# 查看cuffdiff的生成文件。
cd cuffdiff_outls

# 打开 gene_exp.diff文件,按第10列排序,并且只看部分列
cat gene_exp.diff | grep yes | sort -k10n | cut -f 3,8,9,10 | head

Tophat alignment script: run-tophat.sh
[Shell] 纯文本查看 复制代码
# 对每个数据文件运行tophat
# 在错误处停止set -ue
REF=~/refs/hg38/chr4
# 这个文件夹将存放比对结果。
mkdir -p bam
for name in ctl1 ctl2 ctl3 uvb1 uvb2 uvb3
do   
 # 创建fastq文件名。    
FASTQ=data/$name.fastq    
# 输出目录的名字。    
OUTPUT=data/$name-tophat    
# 运行tophat并把结果放到OUTPUT目录    
tophat -G chr4.gtf -o $OUTPUT $REF $FASTQ    
# 将生成的结果文件accepted_hits.bam重新命名。    
mv -f $OUTPUT/accepted_hits.bam bam/$name.bam
done;



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



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






上一篇:Biostar:课程25、26:bedtools的安装和使用
下一篇:Biostar(大结局):课程29、30差异表达分析以及用ErmineJ做GO
回复

使用道具 举报

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

本版积分规则

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

GMT+8, 2019-9-21 14:16 , Processed in 0.035537 second(s), 28 queries .

Powered by Discuz! X3.2

© 2001-2013 Comsenz Inc.