Biostar：课程27、28

公彦栋 2017-10-23

展开全文

翻译：生物女博士

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

# 为作者的注释

#* 为译者的注释

Lecture 27 - RNA-Seq入门

#* 由于课程时间比较上比较旧，所以使用的还是TopHat

#* 目前已经用Hisat2取代了TopHat流程了。具体使用请参考我们公众号文章《新司机带你学RNA-Seq数据分析》by 徐春晖

#* 另外，使用STAR进行RNA-seq分析见《比对软件STAR的使用—高通量测序数据处理学习记录（一）》by 徐春晖

#* 更多RNA-seq分析相关文章可以使用我们公众号全文搜索程序（https://xuzhougeng./biosxy/）（洲更同学荣誉出品）进行搜索。

# 安装TopHat

cd ~/src

# Mac OSX

curl -OL http://ccb./software/tophat/downloads/tophat-2.0.13.OSX_x86_64.tar.gz

# Linux

# curl -OL http://ccb./software/tophat/downloads/tophat-2.0.13.Linux_x86_64.tar.gz

# 解压和安装

tar zxvf tophat-2.0.13.OSX_x86_64.tar.gzln -s ~/src/tophat-2.0.13.OSX_x86_64/tophat ~/bin

# 跟着Stephen Turner的RNA-Seq训练营来学习。

# http://bioconnector./workshops/lessons/rnaseq-1day/#setup

#* 网址已经失效了。找到一个RNAseq分析相关的：http:///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./iua1/courses/files/2014/rnaseq-data.tar.gztar 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.bamsamtools 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./cufflinks/assets/downloads/cufflinks-2.1.1.OSX_x86_64.tar.gztar zxvf cufflinks-2.1.1.OSX_x86_64.tar.gzln -fs ~/src/cufflinks-2.1.1.OSX_x86_64/cufflinks ~/binln -fs ~/src/cufflinks-2.1.1.OSX_x86_64/cuffdiff ~/bin

# Linux下则下载另一个不同版本。注意目录名字的改变。

# 自行根据上边的命令进行修改。

# curl -OL http://cole-trapnell-lab./cufflinks/assets/downloads/cufflinks-2.1.1.Linux_x86_64.tar.gz

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

cuffdiff

# 我们需要用tophat来比对，然后用cuffdiff来进行定量。

# 我们假设基因组已经如lecture27那样，下载并且索引了

#* 这一步获取了一些rnaseq的数据

curl -OL http://www.personal./iua1/courses/files/2014/rnaseq-data.tar.gztar zxvf rnaseq-data.tar.gzrm 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.bamsamtools index bam/ctl1.bam

# 现在来重复之前的每个步骤

# 你应该写一个脚本，或者将所有命令列入一个更自动的脚本（见下文）

# bash run-tophat.sh

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

cuffdiff

# 运行cuffdiff，并把运行结果放到bam/cuffdiff_out文件夹下

cd bamcuffdiff -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

# 对每个数据文件运行tophat# 在错误处停止set -ueREF=~/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;