原创10000+生信教程大神给你的RNA实战视频演练

准备工作

1. 安装conda

推荐使用偷懒方法,比如安装miniconda软件,下载地址:https://mirrors.tuna.tsinghua.edu.cn/anaconda/miniconda/ 这样就可以使用它安装绝大部分其它软件。

但是在中国大陆的小伙伴,需要更改镜像源配置

conda config --add channels https://mirrors.tuna.tsinghua.edu.cn/anaconda/pkgs/free
conda config --add channels https://mirrors.tuna.tsinghua.edu.cn/anaconda/cloud/conda-forge
conda config --add channels https://mirrors.tuna.tsinghua.edu.cn/anaconda/cloud/bioconda
conda config --set show_channel_urls yes

2. 安装软件

为了避免污染linux工作环境,推荐在conda中创建各个流程的安装环境,比如:

conda create -n rna python=2 #创建名为rna的软件安装环境
conda info --envs #查看当前conda环境
source activate rna #激活conda的rna环境

先读文献调研,得到转录组分析需要用到的软件列表;

  • 质控
  • fastqc , multiqc, trimmomatic, cutadapt ,trim-galore
  • 比对
  • star, hisat2, bowtie2, tophat, bwa, subread
  • 计数
  • htseq, bedtools, deeptools, salmon

如果你对一个软件不了解的话,那么安装之前在https://bioconda.github.io/recipes.html,检索该软件包是否存在,或者使用 conda search packagename进行检索。

但是我帮你确定好了下面的软件安装代码是可行的!!!

conda install -y sra-tools
conda install -y trimmomatic
conda install -y cutadapt multiqc 
conda install -y trim-galore
conda install -y star hisat2 bowtie2
conda install -y subread tophat htseq bedtools deeptools
conda install -y salmon
source deactivate #注销当前的rna环境

转录组流程

step1: sra2fastq

下载SRA数据

新建一个名为 SRR_Acc_List.txt的文档,将SRR号码保存在文档内,一个号码占据一行。文件可以在我的GitHub下载获取:https://github.com/jmzeng1314/GEO/blob/master/airway/SRRAccList.txt

  • prefetch下载数据
wkd=/home/jmzeng/project/airway/ #设置工作目录
source activate rna
cat SRR_Acc_List.txt | while read id; do (prefetch  ${id} &);done
ps -ef | grep prefetch | awk '{print $2}' | while read id; do kill ${id}; done #在内地下载速度很慢,所以我杀掉这些下载进程
  • 拷贝事先下载好的sra数据
mkdir $wkd/raw 
cd $wkd/raw 
ls /public/project/RNA/airway/sra/*  | while read id; do ( nohup fastq-dump --gzip --split-3 -O ./ ${id} & ); done
source deactivate
  • 得到的SRA数据如下
/public/project/RNA/airway/sra/
├── [1.6G]  SRR1039508.sra
├── [1.4G]  SRR1039509.sra
├── [1.6G]  SRR1039510.sra
├── [1.5G]  SRR1039511.sra
├── [2.0G]  SRR1039512.sra
├── [2.2G]  SRR1039513.sra
├── [3.0G]  SRR1039514.sra
├── [1.9G]  SRR1039515.sra
├── [2.1G]  SRR1039516.sra
├── [2.6G]  SRR1039517.sra
├── [2.3G]  SRR1039518.sra
├── [2.0G]  SRR1039519.sra
├── [2.1G]  SRR1039520.sra
├── [2.4G]  SRR1039521.sra
├── [2.0G]  SRR1039522.sra
└── [2.2G]  SRR1039523.sra
  • sra格式转fastq格式

格式转还用到的软件是fastq-dump

for i in $wkd/*sra
do
        echo $i
        nohup fastq-dump --split-3 --skip-technical --clip --gzip $i &
done
  • 得到fastq数据如下

原始数据是双端测序结果,fastq-dump配合--split-3参数,一个样本被拆分成两个fastq文件

├── [1.3G]  SRR1039508_1.fastq.gz
├── [1.3G]  SRR1039508_2.fastq.gz
├── [1.2G]  SRR1039509_1.fastq.gz
├── [1.2G]  SRR1039509_2.fastq.gz
├── [1.3G]  SRR1039510_1.fastq.gz
├── [1.3G]  SRR1039510_2.fastq.gz
├── [1.2G]  SRR1039511_1.fastq.gz
├── [1.2G]  SRR1039511_2.fastq.gz
├── [1.6G]  SRR1039512_1.fastq.gz
├── [1.6G]  SRR1039512_2.fastq.gz
├── [950M]  SRR1039513_1.fastq.gz
├── [952M]  SRR1039513_2.fastq.gz
├── [2.4G]  SRR1039514_1.fastq.gz
......
├── [1.5G]  SRR1039522_1.fastq.gz
├── [1.5G]  SRR1039522_2.fastq.gz
├── [1.8G]  SRR1039523_1.fastq.gz
└── [1.8G]  SRR1039523_2.fastq.gz

step2: check quality of sequence reads

fastqc生成质控报告,multiqc将各个样本的质控报告整合为一个。

ls *gz | xargs fastqc -t 10
multiqc ./
  • 得到结果如下
├── [4.0K]  multiqc_data
│   ├── [2.1M]  multiqc_data.json
│   ├── [6.8K]  multiqc_fastqc.txt
│   ├── [2.2K]  multiqc_general_stats.txt
│   ├── [ 16K]  multiqc.log
│   └── [3.4K]  multiqc_sources.txt
├── [1.5M]  multiqc_report.html
├── [236K]  SRR1039508_1_fastqc.html
├── [279K]  SRR1039508_1_fastqc.zip
├── [238K]  SRR1039508_2_fastqc.html
├── [286K]  SRR1039508_2_fastqc.zip
├── [236K]  SRR1039510_1_fastqc.html
├── [278K]  SRR1039510_1_fastqc.zip
├── [241K]  SRR1039510_2_fastqc.html
├── [292K]  SRR1039510_2_fastqc.zip
......
├── [220K]  SRR1039522_fastqc.zip
├── [234K]  SRR1039523_1_fastqc.html
├── [273K]  SRR1039523_1_fastqc.zip
├── [232K]  SRR1039523_2_fastqc.html
└── [274K]  SRR1039523_2_fastqc.zip

每个idfastqc.html都是一个质量报告,multiqcreport.html是所有样本的整合报告

step3: filter the bad quality reads and remove adaptors.

  • 运行如下代码,得到名为config的文件,包含两列数据
mkdir $wkd/clean 
cd $wkd/clean 
ls /home/jmzeng/project/airway/raw/*_1.fastq.gz >1
ls /home/jmzeng/project/airway/raw/*_2.fastq.gz >2
paste 1 2  > config
  • 打开文件 qc.sh ,并且写入如下内容

trim_galore,用于去除低质量和接头数据

source activate rna
bin_trim_galore=trim_galore
dir='/home/jmzeng/project/airway/clean'
cat $1 |while read id
do
        arr=(${id})
        fq1=${arr[0]}
        fq2=${arr[1]} 
nohup $bin_trim_galore -q 25 --phred33 --length 36 --stringency 3 --paired -o $dir  $fq1 $fq2 & 
done 
source deactivate
  • 运行qc.sh
bash qc.sh config #config是传递进去的参数
  • 结果显示如下
├── [2.9K]  SRR1039508_1.fastq.gz_trimming_report.txt
├── [1.2G]  SRR1039508_1_val_1.fq.gz
├── [3.1K]  SRR1039508_2.fastq.gz_trimming_report.txt
├── [1.2G]  SRR1039508_2_val_2.fq.gz
├── [2.9K]  SRR1039509_1.fastq.gz_trimming_report.txt
......
├── [2.9K]  SRR1039522_1.fastq.gz_trimming_report.txt
├── [1.4G]  SRR1039522_1_val_1.fq.gz
├── [3.1K]  SRR1039522_2.fastq.gz_trimming_report.txt
├── [1.4G]  SRR1039522_2_val_2.fq.gz
├── [2.9K]  SRR1039523_1.fastq.gz_trimming_report.txt
├── [1.7G]  SRR1039523_1_val_1.fq.gz
├── [3.1K]  SRR1039523_2.fastq.gz_trimming_report.txt
└── [1.7G]  SRR1039523_2_val_2.fq.gz

step4: alignment

star, hisat2, bowtie2, tophat, bwa, subread都是可以用于比到的软件

  • 先运行一个样本,测试一下
cd $wkd/test 
source activate rna

#比对流程
id=SRR1039508
hisat2 -p 10 -x /public/reference/index/hisat/hg38/genome -1 ${id}1_val_1.fq -2 ${id}2_val_2.fq -S ${id}.hisat.sam
subjunc -T 5 -i /public/reference/index/subread/hg38 -r ${id}1_val_1.fq -R ${id}2_val_2.fq -o ${id}.subjunc.sam  
bowtie2 -p 10 -x /public/reference/index/bowtie/hg38 -1 ${id}1_val_1.fq -2 ${id}2_val_2.fq -S ${id}.bowtie.sam
bwa mem -t 5 -M  /public/reference/index/bwa/hg38 ${id}1_val_1.fq ${id}2_val_2.fq > ${id}.bwa.sam
  • 批量比对代码
cd $wkd/clean 

ls *gz|cut -d"_" -f 1 |sort -u |while read id;do
ls -lh ${id}1_val_1.fq.gz   ${id}2_val_2.fq.gz 
hisat2 -p 10 -x /public/reference/index/hisat/hg38/genome -1 ${id}1_val_1.fq.gz -2 ${id}2_val_2.fq.gz -S ${id}.hisat.sam
subjunc -T 5  -i /public/reference/index/subread/hg38 -r ${id}1_val_1.fq.gz -R ${id}2_val_2.fq.gz -o ${id}.subjunc.sam
bowtie2 -p 10 -x /public/reference/index/bowtie/hg38 -1 ${id}1_val_1.fq.gz -2 ${id}2_val_2.fq.gz -S ${id}.bowtie.sam
bwa mem -t 5 -M /public/reference/index/bwa/hg38 ${id}1_val_1.fq.gz ${id}2_val_2.fq.gz > ${id}.bwa.sam
done
  • sam文件转bam
ls *.sam|while read id ;do (samtools sort -O bam -@ 5  -o $(basename ${id} ".sam").bam   ${id});done
rm *.sam
  • 为bam文件建立索引
ls *.bam |xargs -i samtools index {}
  • reads的比对情况统计
ls *.bam |xargs -i samtools flagstat -@ 10  {}  >
ls *.bam |while read id ;do ( nohup samtools flagstat -@ 1 $id >  $(basename ${id} ".bam").flagstat & );done
source deactivate
  • 最终结果显示如下
├── [1.8G]  SRR1039508.bowite2.bam
├── [2.9M]  SRR1039508.bowite2.bam.bai
├── [ 444]  SRR1039508.bowite2.flagstat
├── [ 10G]  SRR1039508.bowite2.sam
├── [1.7G]  SRR1039509.bowite2.bam
......
├── [2.0G]  SRR1039521.bowite2.bam
├── [2.9M]  SRR1039521.bowite2.bam.bai
├── [ 444]  SRR1039521.bowite2.flagstat
├── [ 10G]  SRR1039521.bowite2.sam
├── [2.3G]  SRR1039522.bowite2.bam
├── [3.0M]  SRR1039522.bowite2.bam.bai
├── [ 444]  SRR1039522.bowite2.flagstat
├── [ 12G]  SRR1039522.bowite2.sam
├── [2.5G]  SRR1039523.bowite2.bam
├── [3.0M]  SRR1039523.bowite2.bam.bai
├── [ 444]  SRR1039523.bowite2.flagstat
└── [ 14G]  SRR1039523.bowite2.sam

step5: counts

mkdir $wkd/align 
cd $wkd/align 
source activate rna
for fn in {508..523}
do
featureCounts -T 5 -p -t exon -g gene_id  -a /public/reference/gtf/gencode/gencode.v25.annotation.gtf.gz -o counts.txt SRR1039$fn.bam
done
source deactivate
  • 得到的文件如下
1 # Program:featureCounts v1.6.1; Command:"featureCounts" "-T" "5" "-p" "-t" "exon" "-g" "gene_id" "-a" "/public/reference/gtf/gencode/ge
      2 Geneid  Chr     Start   End     Strand  Length  /home/llwu/RNA/airway/2.align/bowite2/SRR1039523.bowite2.bam
      3 ENSG00000223972.5       chr1;chr1;chr1;chr1;chr1;chr1;chr1;chr1;chr1    11869;12010;12179;12613;12613;12975;13221;13221;13453   12227;1
      4 ENSG00000227232.5       chr1;chr1;chr1;chr1;chr1;chr1;chr1;chr1;chr1;chr1;chr1  14404;15005;15796;16607;16858;17233;17606;17915;18268;2
      5 ENSG00000278267.1       chr1    17369   17436   -       68      9
      6 ENSG00000243485.4       chr1;chr1;chr1;chr1;chr1;chr1   29554;30267;30366;30564;30976;30976     30039;30667;30503;30667;31097;31109
      7 ENSG00000237613.2       chr1;chr1;chr1;chr1;chr1        34554;35245;35277;35721;35721   35174;35481;35481;36073;36081   -;-;-;-;-
      8 ENSG00000268020.3       chr1    52473   53312   +       840     0
      9 ENSG00000240361.1       chr1    62948   63887   +       940     0
     10 ENSG00000186092.4       chr1    69091   70008   +       918     0

step5: DEG

  • 差异分析之前需要首先对转录组上游分析得到的文件 all.id.txt 进行一定程度的检查,代码如:
rm(list = ls())
options(stringsAsFactors = F)
a=read.table('all.id.txt',header = T)
tmp=a[1:14,1:7]
meta=a[,1:6]
exprSet=a[,7:ncol(a)]
colnames(exprSet)
a2=exprSet[,'SRR1039516.hisat.bam']

library(airway)
data(airway)
exprSet=assay(airway)
colnames(exprSet)
a1=exprSet[,'SRR1039516']
group_list=colData(airway)[,3]

a2=data.frame(id=meta[,1],a2=a2)
a1=data.frame(id=names(a1),a1=as.numeric(a1))
library(stringr)
a2$id <- str_split(a2$id,'\\.',simplify = T)[,1]
tmp=merge(a1,a2,by='id')
png('tmp.png')
plot(tmp[,2:3])
dev.off()

library(corrplot)
png('cor.png')
corrplot(cor(log2(exprSet+1)))
dev.off()

library(pheatmap)
png('heatmap.png')
m=cor(log2(exprSet+1))
pheatmap(scale(cor(log2(exprSet+1))))
dev.off()
  • 检查通过之前就可以正式做差异分析以及后续的富集分析,代码参考:https://github.com/jmzeng1314/my-R/blob/master/10-RNA-seq-3-groups/hisat2mm10htseq.R 其实都在GitHub哦,视频讲解了这个代码如何使用。
  • 更多干货直接扫描看视频吧,还等什么!!!

原文发布于微信公众号 - 生信技能树(biotrainee)

原文发表时间:2018-08-05

本文参与腾讯云自媒体分享计划,欢迎正在阅读的你也加入,一起分享。

发表于

我来说两句

0 条评论
登录 后参与评论

相关文章

来自专栏小文博客

华为Mate10安装黑域

50640
来自专栏逸鹏说道

Stack Overflow 2016最新架构探秘

这篇文章主要揭秘 Stack Overflow 截止到 2016 年的技术架构。   首先给出一个直观的数据,让大家有个初步的印象。   相比于 2013 ...

34070
来自专栏杨建荣的学习笔记

需要了解的pssh(r11笔记第28天)

昨天的一篇文章,关于ssh命令的几个使用小技巧(r11笔记第27天),也收到了不少朋友的反馈,其中有个朋友提议说还是用pssh吧,我想想也是。 对...

35760
来自专栏企鹅号快讯

可编程逻辑控制器PLC

PLC代表可编程逻辑控制器。它们基本上用于控制工业中的自动化系统。它们是最先进和最简单的控制系统之一,现在正在大规模地取代硬连线逻辑继电器。 ? 在详细介绍PL...

26170
来自专栏企鹅号快讯

如何选择TensorFlow安装环境

许多软件都会使用一些库和独立维护的软件包。对于开发者而言,这是一件好事,因为这种做法有利于代码复用,而且他们 可专注于创建新的功能,而无需重复造轮。然而,这种做...

47650
来自专栏linjinhe的专栏

理解 CPU 利用率

在 Linux shell 上执行 top 命令,可以看到这样一行 CPU 利用率的数据:

40060
来自专栏安恒信息

预警 | Windows再爆SMB服务0day漏洞 安恒信息提供一键修复工具

在刚刚结束的2017年度DEFCON大会上,安全研究人员曝光了一个潜伏在Windows系统上长达20年之久的安全漏洞。该漏洞被命名为SMBLoris,攻击者使用...

362110
来自专栏杨海春的专栏

常用机器性能评估工具

软件系统跑在机器上,处理能力受硬件制约,所以,单机处理能力会有上限。评估机器处理能力的上限,检查程序的瓶颈在哪,有助于程序性能分析。主要的几大硬件:CPU、内存...

69000
来自专栏数据和云

ASM的备份解析与恢复

编辑说明:《Oracle性能优化与诊断案例精选》出版以来,收到很多读者的来信和评论,我们会通过连载的形式将书中内容公布出来,希望书中内容能够帮助到更多的读者朋友...

34160
来自专栏IT技术精选文摘

缓存更新的套路

看到好些人在写更新缓存数据代码时,先删除缓存,然后再更新数据库,而后续的操作会把数据再装载的缓存中。然而,这个是逻辑是错误的。试想,两个并发操作,一个是更新操作...

44070

扫码关注云+社区

领取腾讯云代金券