Biostar：课程27、28
翻译：生物女博士
注：本系列课程翻译已获得原作者授权。
# 为作者的注释
#* 为译者的注释
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
curl -OL http://ccb.jhu.edu/software/tophat/downloads/tophat-2.0.13.OSX_x86_64.tar.gz# Linux
# curl -OL http://ccb.jhu.edu/software/tophat/downloads/tophat-2.0.13.Linux_x86_64.tar.gz# 解压和安装
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训练营来学习。
# http://bioconnector.github.io/workshops/lessons/rnaseq-1day/#setup
#* 网址已经失效了。找到一个RNAseq分析相关的：http://bioconnector.org/workshops/r-rnaseq-airway.html
# 从Ensembl下载基因组数据。
# 为人类基因组创建一个目录
mkdir -p ~/refs/hg38# 获取基因组序列
curl -L ftp://ftp.ensembl.org/pub/release-77/fasta/homo_sapiens/dna/Homo_sapiens.GRCh38.dna.chromosome.4.fa.gz > ~/refs/hg38/chr4.fa.gz# 获取注释信息
curl -L ftp://ftp.ensembl.org/pub/release-77/gtf/homo_sapiens/Homo_sapiens.GRCh38.77.gtf.gz > ~/refs/hg38/hg38.gtf.gz# 解压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# 获取和解压数据。
curl -OL http://www.personal.psu.edu/iua1/courses/files/2014/rnaseq-data.tar.gz
tar zxvf rnaseq-data.tar.gz# 首先，我们来比对这个数据
# 给从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.io/cufflinks/assets/downloads/cufflinks-2.1.1.OSX_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下则下载另一个不同版本。注意目录名字的改变。
# 自行根据上边的命令进行修改。
# curl -OL http://cole-trapnell-lab.github.io/cufflinks/assets/downloads/cufflinks-2.1.1.Linux_x86_64.tar.gz
# 来看一下是不是可以正常使用了。
cuffdiff# 我们需要用tophat来比对，然后用cuffdiff来进行定量。
# 我们假设基因组已经如lecture27那样，下载并且索引了
#* 这一步获取了一些rnaseq的数据
curl -OL http://www.personal.psu.edu/iua1/courses/files/2014/rnaseq-data.tar.gz
tar zxvf rnaseq-data.tar.gz
rm rnaseq-data.tar.gz# 如果没有，那就通过下边链接获取。
# curl -L ftp://ftp.ensembl.org/pub/release-77/gtf/homo_sapiens/Homo_sapiens.GRCh38.77.gtf.gz > ~/refs/hg38/hg38.gtf.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_out
ls# 打开 gene_exp.diff文件，按第10列排序，并且只看部分列
cat gene_exp.diff | grep yes | sort -k10n | cut -f 3,8,9,10 | headTophat alignment script: run-tophat.sh
# 对每个数据文件运行tophat# 在错误处停止set -ue
REF=~/refs/hg38/chr4# 这个文件夹将存放比对结果。mkdir -p bamfor name in ctl1 ctl2 ctl3 uvb1 uvb2 uvb3do
# 创建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.bamdone;Biostar非常适合有生物学基础且想自学生信的这部分人入门。译者正是通过该课程入门的生信。
可点击 阅读原文 或在公众号的 菜单→课程系列 中找到本系列课程。
欢迎关注我们