🔥 热搜 🔥
上海习近平新疆鄂州父女瓜乌鲁木齐疫情H工口小学生赛高习明泽芊川一笑图包印尼排华
分类
社会娱乐国际人权科技经济其它
🔥 热搜 🔥
上海习近平新疆鄂州父女瓜乌鲁木齐疫情H工口小学生赛高习明泽芊川一笑图包印尼排华
分类
社会娱乐国际人权科技经济其它
查看原文
其他

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

生信技能树 生信技能树 2022-06-07
收录于合集 #RNA 36个

准备工作

1. 安装conda

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

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

  1. conda config --add channels https://mirrors.tuna.tsinghua.edu.cn/anaconda/pkgs/free

  2. conda config --add channels https://mirrors.tuna.tsinghua.edu.cn/anaconda/cloud/conda-forge

  3. conda config --add channels https://mirrors.tuna.tsinghua.edu.cn/anaconda/cloud/bioconda

  4. conda config --set show_channel_urls yes

2. 安装软件

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

  1. conda create -n rna python=2 #创建名为rna的软件安装环境

  2. conda info --envs #查看当前conda环境

  3. 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进行检索。

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

  1. conda install -y sra-tools

  2. conda install -y trimmomatic

  3. conda install -y cutadapt multiqc

  4. conda install -y trim-galore

  5. conda install -y star hisat2 bowtie2

  6. conda install -y subread tophat htseq bedtools deeptools

  7. conda install -y salmon

  8. source deactivate #注销当前的rna环境

转录组流程

step1: sra2fastq

下载SRA数据

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

  • prefetch下载数据

  1. wkd=/home/jmzeng/project/airway/ #设置工作目录

  2. source activate rna

  3. cat SRR_Acc_List.txt | while read id; do (prefetch  ${id} &);done

  4. ps -ef | grep prefetch | awk '{print $2}' | while read id; do kill ${id}; done #在内地下载速度很慢，所以我杀掉这些下载进程

  • 拷贝事先下载好的sra数据

  1. mkdir $wkd/raw

  2. cd $wkd/raw

  3. ls /public/project/RNA/airway/sra/*  | while read id; do ( nohup fastq-dump --gzip --split-3 -O ./ ${id} & ); done

  4. source deactivate

  • 得到的SRA数据如下

  1. /public/project/RNA/airway/sra/

  2. ├── [1.6G]  SRR1039508.sra

  3. ├── [1.4G]  SRR1039509.sra

  4. ├── [1.6G]  SRR1039510.sra

  5. ├── [1.5G]  SRR1039511.sra

  6. ├── [2.0G]  SRR1039512.sra

  7. ├── [2.2G]  SRR1039513.sra

  8. ├── [3.0G]  SRR1039514.sra

  9. ├── [1.9G]  SRR1039515.sra

  10. ├── [2.1G]  SRR1039516.sra

  11. ├── [2.6G]  SRR1039517.sra

  12. ├── [2.3G]  SRR1039518.sra

  13. ├── [2.0G]  SRR1039519.sra

  14. ├── [2.1G]  SRR1039520.sra

  15. ├── [2.4G]  SRR1039521.sra

  16. ├── [2.0G]  SRR1039522.sra

  17. └── [2.2G]  SRR1039523.sra

  • sra格式转fastq格式

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

  1. for i in $wkd/*sra

  2. do

  3.        echo $i

  4.        nohup fastq-dump --split-3 --skip-technical --clip --gzip $i &

  5. done

  • 得到fastq数据如下

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

  1. ├── [1.3G]  SRR1039508_1.fastq.gz

  2. ├── [1.3G]  SRR1039508_2.fastq.gz

  3. ├── [1.2G]  SRR1039509_1.fastq.gz

  4. ├── [1.2G]  SRR1039509_2.fastq.gz

  5. ├── [1.3G]  SRR1039510_1.fastq.gz

  6. ├── [1.3G]  SRR1039510_2.fastq.gz

  7. ├── [1.2G]  SRR1039511_1.fastq.gz

  8. ├── [1.2G]  SRR1039511_2.fastq.gz

  9. ├── [1.6G]  SRR1039512_1.fastq.gz

  10. ├── [1.6G]  SRR1039512_2.fastq.gz

  11. ├── [950M]  SRR1039513_1.fastq.gz

  12. ├── [952M]  SRR1039513_2.fastq.gz

  13. ├── [2.4G]  SRR1039514_1.fastq.gz

  14. ......

  15. ├── [1.5G]  SRR1039522_1.fastq.gz

  16. ├── [1.5G]  SRR1039522_2.fastq.gz

  17. ├── [1.8G]  SRR1039523_1.fastq.gz

  18. └── [1.8G]  SRR1039523_2.fastq.gz

step2: check quality of sequence reads

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

  1. ls *gz | xargs fastqc -t 10

  2. multiqc ./

  • 得到结果如下

  1. ├── [4.0K]  multiqc_data

  2. │   ├── [2.1M]  multiqc_data.json

  3. │   ├── [6.8K]  multiqc_fastqc.txt

  4. │   ├── [2.2K]  multiqc_general_stats.txt

  5. │   ├── [ 16K]  multiqc.log

  6. │   └── [3.4K]  multiqc_sources.txt

  7. ├── [1.5M]  multiqc_report.html

  8. ├── [236K]  SRR1039508_1_fastqc.html

  9. ├── [279K]  SRR1039508_1_fastqc.zip

  10. ├── [238K]  SRR1039508_2_fastqc.html

  11. ├── [286K]  SRR1039508_2_fastqc.zip

  12. ├── [236K]  SRR1039510_1_fastqc.html

  13. ├── [278K]  SRR1039510_1_fastqc.zip

  14. ├── [241K]  SRR1039510_2_fastqc.html

  15. ├── [292K]  SRR1039510_2_fastqc.zip

  16. ......

  17. ├── [220K]  SRR1039522_fastqc.zip

  18. ├── [234K]  SRR1039523_1_fastqc.html

  19. ├── [273K]  SRR1039523_1_fastqc.zip

  20. ├── [232K]  SRR1039523_2_fastqc.html

  21. └── [274K]  SRR1039523_2_fastqc.zip

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

step3: filter the bad quality reads and remove adaptors.

  • 运行如下代码，得到名为config的文件，包含两列数据

  1. mkdir $wkd/clean

  2. cd $wkd/clean

  3. ls /home/jmzeng/project/airway/raw/*_1.fastq.gz >1

  4. ls /home/jmzeng/project/airway/raw/*_2.fastq.gz >2

  5. paste 1 2  > config

  • 打开文件 qc.sh ，并且写入如下内容

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

  1. source activate rna

  2. bin_trim_galore=trim_galore

  3. dir='/home/jmzeng/project/airway/clean'

  4. cat $1 |while read id

  5. do

  6.        arr=(${id})

  7.        fq1=${arr[0]}

  8.        fq2=${arr[1]}

  9. nohup $bin_trim_galore -q 25 --phred33 --length 36 --stringency 3 --paired -o $dir  $fq1 $fq2 &

  10. done

  11. source deactivate

  • 运行qc.sh

  1. bash qc.sh config #config是传递进去的参数

  • 结果显示如下

  1. ├── [2.9K]  SRR1039508_1.fastq.gz_trimming_report.txt

  2. ├── [1.2G]  SRR1039508_1_val_1.fq.gz

  3. ├── [3.1K]  SRR1039508_2.fastq.gz_trimming_report.txt

  4. ├── [1.2G]  SRR1039508_2_val_2.fq.gz

  5. ├── [2.9K]  SRR1039509_1.fastq.gz_trimming_report.txt

  6. ......

  7. ├── [2.9K]  SRR1039522_1.fastq.gz_trimming_report.txt

  8. ├── [1.4G]  SRR1039522_1_val_1.fq.gz

  9. ├── [3.1K]  SRR1039522_2.fastq.gz_trimming_report.txt

  10. ├── [1.4G]  SRR1039522_2_val_2.fq.gz

  11. ├── [2.9K]  SRR1039523_1.fastq.gz_trimming_report.txt

  12. ├── [1.7G]  SRR1039523_1_val_1.fq.gz

  13. ├── [3.1K]  SRR1039523_2.fastq.gz_trimming_report.txt

  14. └── [1.7G]  SRR1039523_2_val_2.fq.gz

step4: alignment

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

  • 先运行一个样本，测试一下

  1. cd $wkd/test

  2. source activate rna

  4. #比对流程

  5. id=SRR1039508

  6. 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

  7. 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  

  8. 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

  9. bwa mem -t 5 -M  /public/reference/index/bwa/hg38 ${id}1_val_1.fq ${id}2_val_2.fq > ${id}.bwa.sam

  • 批量比对代码

  1. cd $wkd/clean

  3. ls *gz|cut -d"_" -f 1 |sort -u |while read id;do

  4. ls -lh ${id}1_val_1.fq.gz   ${id}2_val_2.fq.gz

  5. 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

  6. 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

  7. 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

  8. 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

  9. done

  • sam文件转bam

  1. ls *.sam|while read id ;do (samtools sort -O bam -@ 5  -o $(basename ${id} ".sam").bam   ${id});done

  2. rm *.sam

  • 为bam文件建立索引

  1. ls *.bam |xargs -i samtools index {}

  • reads的比对情况统计

  1. ls *.bam |xargs -i samtools flagstat -@ 10  {}  >

  2. ls *.bam |while read id ;do ( nohup samtools flagstat -@ 1 $id >  $(basename ${id} ".bam").flagstat & );done

  3. source deactivate

  • 最终结果显示如下

  1. ├── [1.8G]  SRR1039508.bowite2.bam

  2. ├── [2.9M]  SRR1039508.bowite2.bam.bai

  3. ├── [ 444]  SRR1039508.bowite2.flagstat

  4. ├── [ 10G]  SRR1039508.bowite2.sam

  5. ├── [1.7G]  SRR1039509.bowite2.bam

  6. ......

  7. ├── [2.0G]  SRR1039521.bowite2.bam

  8. ├── [2.9M]  SRR1039521.bowite2.bam.bai

  9. ├── [ 444]  SRR1039521.bowite2.flagstat

  10. ├── [ 10G]  SRR1039521.bowite2.sam

  11. ├── [2.3G]  SRR1039522.bowite2.bam

  12. ├── [3.0M]  SRR1039522.bowite2.bam.bai

  13. ├── [ 444]  SRR1039522.bowite2.flagstat

  14. ├── [ 12G]  SRR1039522.bowite2.sam

  15. ├── [2.5G]  SRR1039523.bowite2.bam

  16. ├── [3.0M]  SRR1039523.bowite2.bam.bai

  17. ├── [ 444]  SRR1039523.bowite2.flagstat

  18. └── [ 14G]  SRR1039523.bowite2.sam

step5: counts

  1. mkdir $wkd/align

  2. cd $wkd/align

  3. source activate rna

  4. for fn in {508..523}

  5. do

  6. 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

  7. done

  8. source deactivate

  • 得到的文件如下

  1.      1 # Program:featureCounts v1.6.1; Command:"featureCounts" "-T" "5" "-p" "-t" "exon" "-g" "gene_id" "-a" "/public/reference/gtf/gencode/ge

  2.      2 Geneid  Chr     Start   End     Strand  Length  /home/llwu/RNA/airway/2.align/bowite2/SRR1039523.bowite2.bam

  3.      3 ENSG00000223972.5       chr1;chr1;chr1;chr1;chr1;chr1;chr1;chr1;chr1    11869;12010;12179;12613;12613;12975;13221;13221;13453   12227;1

  4.      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.      5 ENSG00000278267.1       chr1    17369   17436   -       68      9

  6.      6 ENSG00000243485.4       chr1;chr1;chr1;chr1;chr1;chr1   29554;30267;30366;30564;30976;30976     30039;30667;30503;30667;31097;31109

  7.      7 ENSG00000237613.2       chr1;chr1;chr1;chr1;chr1        34554;35245;35277;35721;35721   35174;35481;35481;36073;36081   -;-;-;-;-

  8.      8 ENSG00000268020.3       chr1    52473   53312   +       840     0

  9.      9 ENSG00000240361.1       chr1    62948   63887   +       940     0

  10.     10 ENSG00000186092.4       chr1    69091   70008   +       918     0

step5: DEG

  • 差异分析之前需要首先对转录组上游分析得到的文件 all.id.txt 进行一定程度的检查，代码如：

  1. rm(list = ls())

  2. options(stringsAsFactors = F)

  3. a=read.table('all.id.txt',header = T)

  4. tmp=a[1:14,1:7]

  5. meta=a[,1:6]

  6. exprSet=a[,7:ncol(a)]

  7. colnames(exprSet)

  8. a2=exprSet[,'SRR1039516.hisat.bam']

  10. library(airway)

  11. data(airway)

  12. exprSet=assay(airway)

  13. colnames(exprSet)

  14. a1=exprSet[,'SRR1039516']

  15. group_list=colData(airway)[,3]

  17. a2=data.frame(id=meta[,1],a2=a2)

  18. a1=data.frame(id=names(a1),a1=as.numeric(a1))

  19. library(stringr)

  20. a2$id <- str_split(a2$id,'\\.',simplify = T)[,1]

  21. tmp=merge(a1,a2,by='id')

  22. png('tmp.png')

  23. plot(tmp[,2:3])

  24. dev.off()

  26. library(corrplot)

  27. png('cor.png')

  28. corrplot(cor(log2(exprSet+1)))

  29. dev.off()

  31. library(pheatmap)

  32. png('heatmap.png')

  33. m=cor(log2(exprSet+1))

  34. pheatmap(scale(cor(log2(exprSet+1))))

  35. dev.off()

  • 检查通过之前就可以正式做差异分析以及后续的富集分析，代码参考：https://github.com/jmzeng1314/my-R/blob/master/10-RNA-seq-3-groups/hisat2mm10htseq.R 其实都在GitHub哦，视频讲解了这个代码如何使用。

  • 更多干货直接扫描看视频吧，还等什么！！！

如果需要组装自己的服务器；代办生物信息学服务器

如果需要帮忙下载海外数据(GEO/TCGA/GTEx等等)，点我？

如果需要线下辅导及培训，看招学徒

如果需要个人电脑：个人计算机推荐

如果需要置办生物信息学书籍，看：生信人必备书单

如果需要实习岗位：实习职位发布

价值999的全外显子数据实战基础教程限时免费领取： 外显子福利

价值599的GEO数据库挖掘手把手教程限时免费领取：GEO挖掘福利

更多干货代码，点击下面的阅读原文直达哦！

您可能也对以下帖子感兴趣

文章有问题？点此查看未经处理的缓存