RNA-seq

数据来源https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE49378
相关文献链接https://bmcgenomics.biomedcentral.com/articles/10.1186/1471-2164-14-578#Sec2
通过用RNA-Seq比较自养育空地区、山东以及加拿大育空地区的植物样本来比较不同环境下Eutrema的基因表达模式，文献出用到了9个样本，这里我取其中的两个进行具体分析不同环境下基因的差异表达及相关信息。
样本一：来自加拿大育空地区Yukon Field
SRR1004972
样本二：来自山东Shandong Cabinet
SRR1004966
参考基因组来自Eutrema salsugineum

1、下载数据

mkdir RNAseq_Eutrema_salsugineum #建立一个文件夹
cd RNAseq_Eutrema_salsugineum
prefetch SRR1004972
#出现了error，没有成功，然后尝试用aspera下载
~/.aspera/connect/bin/ascp -T -i /home/yuan/.aspera/connect/etc/asperaweb_id_dsa.openssh -k 1 -l 200m anonftp@ftp-private.ncbi.nlm.nih.gov:/sra/sra-instant/reads/ByRun/sra/SRR/SRR100/SRR1004972/SRR1004972.sra ./
#还是出现了问题
#最终尝试在EMBL-ENA（https://www.ebi.ac.uk/ena/browser/）里通过SRA号找到下载地址，然后用wget下载文件，可以下载，不过下载速度特别慢
wget ftp://ftp.sra.ebi.ac.uk/vol1/fastq/SRR100/006/SRR1004966/SRR1004966.fastq.gz
wget ftp://ftp.sra.ebi.ac.uk/vol1/fastq/SRR100/002/SRR1004972/SRR1004972.fastq.gz
#参考基因组文件
 wget ftp://ftp.ncbi.nlm.nih.gov/genomes/all/GCA/000/478/725/GCA_000478725.1_Eutsalg1_0/GCA_000478725.1_Eutsalg1_0_rna_from_genomic.fna.gz
#解压
gzip -d GCA_000478725.1_Eutsalg1_0_rna_from_genomic.fna.gz
#注释文件
 wget ftp://ftp.ncbi.nlm.nih.gov/genomes/all/GCA/000/478/725/GCA_000478725.1_Eutsalg1_0/GCA_000478725.1_Eutsalg1_0_genomic.gff.gz
gzip -d GCA_000478725.1_Eutsalg1_0_genomic.gff.gz

image.png

2、原始数据质控

fastqc SRR1004966.fastq.gz
fastqc SRR1004972.fastq.gz

3、bwa建立索引

bwa index -a bwtsw GCA_000478725.1_Eutsalg1_0_rna_from_genomic.fna

image.png

4、利用比对软件 bwa 将短序列比对到参考基因组

bwa mem -t 2 -M GCA_000478725.1_Eutsalg1_0_rna_from_genomic.fna SRR1004966.fastq.gz > SRR1004966.sam
bwa mem -t 2 -M GCA_000478725.1_Eutsalg1_0_rna_from_genomic.fna SRR1004972.fastq.gz > SRR1004972.sam

image.png

5、利用 samtools 对 sam 格式的比对文件进行处理，以便进行后续分析

(1)将 sam 格式转化成二进制格式 bam

 samtools view -bt GCA_000478725.1_Eutsalg1_0_rna_from_genomic.fna -@ 2 -o SRR1004966.bam SRR1004966.sam 2>>samtools.log
 samtools view -bt GCA_000478725.1_Eutsalg1_0_rna_from_genomic.fna -@ 2 -o SRR1004972.bam SRR1004972.sam 2>>samtools.log

(2)对 bam 文件中的内容进行排序

samtools sort -O bam -@ 2 -o SRR1004966.sort.bam -T tmp_samtools SRR1004966.bam 2>>samtools.log
samtools sort -O bam -@ 2 -o SRR1004972.sort.bam -T tmp_samtools SRR1004972.bam 2>>samtools.log

(3)为排序后的bam文件建立索引

samtools index SRR1004966.sort.bam
samtools index SRR1004972.sort.bam

(4)显示基因组比对情况

 samtools tview SRR1004966.sort.bam GCA_000478725.1_Eutsalg1_0_rna_from_genomic.fna
 samtools tview SRR1004972.sort.bam GCA_000478725.1_Eutsalg1_0_rna_from_genomic.fna

image.png

(5)测试参考基因组每个位点或一段区域的测序深度

samtools depth SRR1004972.sort.bam >>depth.txt
samtools depth SRR1004966.sort.bam >>depth.txt
less depth.txt

image.png

(6)统计比对结果

samtools flagstat SRR1004972.sort.bam
samtools flagstat SRR1004966.sort.bam

image.png

6、利用 htsep-count进行基因表达量计算， 计算比对到每个基因的短序列数目

htseq-count -f sam -s no -t CDS -i ID -m union -o SRR1004966_assigned.sam SRR1004966.sam GCA_000478725.1_Eutsalg1_0_genomic.gff > SRR1004966_expression.counts
htseq-count -f sam -s no -t CDS -i ID -m union -o SRR1004972_assigned.sam SRR1004972.sam GCA_000478725.1_Eutsalg1_0_genomic.gff > SRR1004972_expression.counts

image.png

7、用DESeq2进行基因表达差异分析

(1)安装R

因为之前在课堂上在windows系统上安装过R，先尝试一下能否安装R包。

image.png

出现了一些问题，又下载并安装了Rtools，问题还是没有解决。
尝试在Linux上安装，经历很多次失败之后终于装上了R version 3.4.4

sudo apt-key adv --keyserver keyserver.ubuntu.com --recv-keys A507B2BBA7803E3B    #A507B2BBA7803E3B为公钥
sudo apt-get update
sudo apt-get install r-base

image.png

(2)执行命令R 进入R环境,并安装差异表达分析包 DESeq2

R
>install.packages("BiocManager")
>BiocManager::install('DEseq2')
>BiocManager::install('ggplots') #画图工具包
>BiocManager::install('clusterProfiler'） #KEGG、GO富集分析工具包

又出现了问题，尝试了很多种方法，一直安装不成功，实验被迫终止，后期问题如果能够解决，会再更新。

尝试写了一下脚本，由于能力和知识不够，所以只针对本实验，相当于一个流程的总结。
执行前应安装好一下软件：
Aspera3.6.2
java 1.8.0_172
FastQC v0.11.7
samtools 1.7
BWA 0.7.17-r1188
HTSeq-0.9.1

#!/bin/bash
mkdir RNAseq_Eutrema_salsugineum  #建立一个文件夹
cd RNAseq_Eutrema_salsugineum
#在EMBL-ENA（https://www.ebi.ac.uk/ena/browser/）里通过SRA号找到下载地址，然后用wget下载文件，可以下载，不过下载速度特别慢
wget ftp://ftp.sra.ebi.ac.uk/vol1/fastq/SRR100/006/SRR1004966/SRR1004966.fastq.gz
wget ftp://ftp.sra.ebi.ac.uk/vol1/fastq/SRR100/002/SRR1004972/SRR1004972.fastq.gz
wget ftp://ftp.ncbi.nlm.nih.gov/genomes/all/GCA/000/478/725/GCA_000478725.1_Eutsalg1_0/GCA_000478725.1_Eutsalg1_0_rna_from_genomic.fna.gz#下载参考基因组文件
gzip -d GCA_000478725.1_Eutsalg1_0_rna_from_genomic.fna.gz#解压
 wget ftp://ftp.ncbi.nlm.nih.gov/genomes/all/GCA/000/478/725/GCA_000478725.1_Eutsalg1_0/GCA_000478725.1_Eutsalg1_0_genomic.gff.gz
gzip -d GCA_000478725.1_Eutsalg1_0_genomic.gff.gz#下载注释文件并解压
fastqc SRR1004966.fastq.gz
fastqc SRR1004972.fastq.gz#原始数据质控
bwa index -a bwtsw GCA_000478725.1_Eutsalg1_0_rna_from_genomic.fna#bwa建立索引
bwa mem -t 2 -M GCA_000478725.1_Eutsalg1_0_rna_from_genomic.fna SRR1004966.fastq.gz > SRR1004966.sam
bwa mem -t 2 -M GCA_000478725.1_Eutsalg1_0_rna_from_genomic.fna SRR1004972.fastq.gz > SRR1004972.sam#利用比对软件 bwa 将短序列比对到参考基因组
#利用 samtools 对 sam 格式的比对文件进行处理，以便进行后续分析
samtools view -bt GCA_000478725.1_Eutsalg1_0_rna_from_genomic.fna -@ 2 -o SRR1004966.bam SRR1004966.sam 2>>samtools.log
 samtools view -bt GCA_000478725.1_Eutsalg1_0_rna_from_genomic.fna -@ 2 -o SRR1004972.bam SRR1004972.sam 2>>samtools.log#将 sam 格式转化成二进制格式 bam
samtools sort -O bam -@ 2 -o SRR1004966.sort.bam -T tmp_samtools SRR1004966.bam 2>>samtools.log
samtools sort -O bam -@ 2 -o SRR1004972.sort.bam -T tmp_samtools SRR1004972.bam 2>>samtools.log#对 bam 文件中的内容进行排序
samtools index SRR1004966.sort.bam
samtools index SRR1004972.sort.bam#为排序后的bam文件建立索引
samtools tview SRR1004966.sort.bam GCA_000478725.1_Eutsalg1_0_rna_from_genomic.fna
samtools tview SRR1004972.sort.bam GCA_000478725.1_Eutsalg1_0_rna_from_genomic.fna#显示基因组比对情况
samtools depth SRR1004972.sort.bam >>depth.txt
samtools depth SRR1004966.sort.bam >>depth.txt
less depth.txt#测试参考基因组每个位点或一段区域的测序深度
samtools flagstat SRR1004972.sort.bam
samtools flagstat SRR1004966.sort.bam#统计比对结果
htseq-count -f sam -s no -t CDS -i ID -m union -o SRR1004966_assigned.sam SRR1004966.sam GCA_000478725.1_Eutsalg1_0_genomic.gff > SRR1004966_expression.counts
htseq-count -f sam -s no -t CDS -i ID -m union -o SRR1004972_assigned.sam SRR1004972.sam GCA_000478725.1_Eutsalg1_0_genomic.gff > SRR1004972_expression.counts#利用 htsep-count进行基因表达量计算， 计算比对到每个基因的短序列数目
#后面得到的数据就可以用来进入R环境，用DESeq2进行基因表达差异分析。