RNA-seq流程学习笔记（9）-使用RSeQC软件对生成的BAM文件进行质控

2023-11-06 11:10
文章标签 进行 流程 学习 使用 笔记 软件 生成 seq bam rna 质控 rseqc

本文主要是介绍RNA-seq流程学习笔记（9）-使用RSeQC软件对生成的BAM文件进行质控，希望对大家解决编程问题提供一定的参考价值，需要的开发者们随着小编来一起学习吧！

参考文章：
用RSeQC对比对后的转录组数据进行质控
高通量测序质控及可视化工具包RSeQC
RSeQC使用笔记

1. 质控的原因及相关软件

在A survey of best practices for RNA-seq data analysis里面，提到了人类基因组应该有70%~90%的比对率，并且多比对read（multi-mapping reads）数量要少。另外比对在外显子和所比对链（uniformity of read coverage on exons and the mapped strand）的覆盖度要保持一致。因此，可以对之前得到的BAM比对文件进行质检。

对BAM文件进行QC的软件包括：

Qualimap：对二代数据进行质控的综合软件

Picard：综合质控学习软件。

RSeQC是发表于2012年的一个RNA-Seq质控工具，属于python包。提供了一系列有用的小工具能够评估高通量测序尤其是RNA-seq数据。比如一些基本模块：检查序列质量、核酸组分偏性、PCR偏性、GC含量偏性，还有RNA-seq特异性模块：评估测序饱和度、映射读数分布、覆盖均匀性、链特异性、转录水平RNA完整性等。

2. RSeQC软件安装

参照文章：RNA-seq流程学习笔记（3）
查看Conda官网Index的RSeQC软件介绍，发现支持python3.6版本，因此直接使用Miniconda3安装， 安装完成后并没有RSeQC这个软件，而是增加了一些python命令，如下：
在这里插入图片描述
虽然该软件的使用命令非常多，但很多功能并不是用来诊断转录组测序的，所以不在我们的考虑范围内。

3.RSeQC处理4种文件格式:

  1. BED 格式：Tab 分割，12列的表示基因模型的纯文本文件。
  2. SAM 或BAM 格式： 用来存储reads比对结果信息，SAM是可读的纯文本文件，然而BAM是SAM的二进制文本，一个压缩的可索引的reads比对文件。
  3. 染色体大小文件： 只有两列的纯文本文件，在“生物信息学文本处理大杂烩(一)”里已经讲过。hg19.chrom_24.sizes是人基因组hg19版本的size文件，是使用UCSC 的fetchChromSizes下载的。
  4. Fasta文件。
    我主要使用的是比对后得到的BAM格式文件。

4. RSeQC软件进行质控检测

1. 使用bam_stat.py命令查看比对的总体情况

#命令说明
Usage: bam_stat.py [options]
Summarizing mapping statistics of a BAM or SAM file.
Options:--version             show program's version number and exit-h, --help            show this help message and exit-i INPUT_FILE, --input-file=INPUT_FILEAlignment file in BAM or SAM format.-q MAP_QUAL, --mapq=MAP_QUALMinimum mapping quality (phred scaled) to determine"uniquely mapped" reads. default=30

#操作记录
(base) zexing@DNA:~/projects/zhaoxiujuan/aligned$ bam_stat.py -i Scr.bam.sort
Load BAM file ...  Done#==================================================
#All numbers are READ count
#==================================================Total records:                          50976263QC failed:                              0
Optical/PCR duplicate:                  0
Non primary hits                        5208051   #表示多匹配位点
Unmapped reads:                         1232377
mapq < mapq_cut (non-unique):           2464685
mapq >= mapq_cut (unique):              42071150
Read-1:                                 21096560
Read-2:                                 20974590
Reads map to '+':                       21026875
Reads map to '-':                       21044275
Non-splice reads:                       18498057
Splice reads:                           23573093
Reads mapped in proper pairs:           41139130
Proper-paired reads map to different chrom:0

确认代码：

bam_stat.py -i ${i}.bam.sort

2. 使用read_distribution.py命令查看基因组覆盖率

#命令说明：
Usage: read_distribution.py [options]
Check reads distribution over exon, intron, UTR, intergenic ... etc
The following reads will be skipped:qc_failedPCR duplicateUnmappedNon-primary (or secondary)
Options:--version             show program's version number and exit-h, --help            show this help message and exit-i INPUT_FILE, --input-file=INPUT_FILEAlignment file in BAM or SAM format.-r REF_GENE_MODEL, --refgene=REF_GENE_MODELReference gene model in bed format.
#该命令需要输入两个文件
# -i为BAM或SAM文件
# -r为参考的bed文件

BED文件可以直接从RSeQC网站下载：
此次练习中下载了人的BED文件，测序数据为小鼠来源，错误！！
在这里插入图片描述下载如下：

#使用wget下载，发现报错，暂未解决。
#在windows 10下面右键保存（2.6M）
#使用scp file name@222.22.222.222：~/dir 命令将文件上传至服务器
#使用gzip命令对文件解压缩，压缩文件不能被识别
(base) zexing@DNA:~/projects/zhaoxiujuan/aligned/BED_file$ gzip -d hg19_RefSeq.bed.gz
(base) zexing@DNA:~/projects/zhaoxiujuan/aligned/BED_file$ ls
hg19_RefSeq.bed

操作记录：

(base) zexing@DNA:~/projects/zhaoxiujuan/aligned$ read_distribution.py -i ./Scr.bam.sort -r /f/xudonglab/zexing/projects/zhaoxiujuan/aligned/BED_file/hg19_RefSeq.bed
processing /f/xudonglab/zexing/projects/zhaoxiujuan/aligned/BED_file/hg19_RefSeq.bed ... Done
processing ./Scr.bam.sort ... FinishedTotal Reads                   44535835
Total Tags                    75880491
Total Assigned Tags           37717256
=====================================================================
Group               Total_bases         Tag_count           Tags/Kb
CDS_Exons           35271889            808827              22.93
5'UTR_Exons         13156148            236757              18.00
3'UTR_Exons         35031450            586598              16.74
Introns             1249039112          30136632            24.13
TSS_up_1kb          19566867            461946              23.61
TSS_up_5kb          88704854            1929741             21.75
TSS_up_10kb         163160166           3349520             20.53
TES_down_1kb        20819429            527166              25.32
TES_down_5kb        89844112            1590630             17.70
TES_down_10kb       160449429           2598922             16.20
=====================================================================

可以用一个饼图来表示，在生信技能树论坛里面还有人专门提问过。（有待进一步研究）
关于RSeQC输出结果的保存，可以使用定向写入【>】来保存。

#使用重定向>命令将read_distribution命令的输出结果定向到指定Log文件中
(base) zexing@DNA:~/projects/zhaoxiujuan/aligned$ read_distribution.py -i ./Scr.bam.sort -r /f/xudonglab/zexing/projects/zhaoxiujuan/aligned/BED_file/hg19_RefSeq.bed > Scr_distribution.log
processing /f/xudonglab/zexing/projects/zhaoxiujuan/aligned/BED_file/hg19_RefSeq.bed ... Done
processing ./Scr.bam.sort ... Finished
(base) zexing@DNA:~/projects/zhaoxiujuan/aligned$ ll
total 277G
-rw-rw-r-- 1 zexing zexing 1.2K 6月   4 13:05 Scr_distribution.log
#使用cat命令查看运行结果
(base) zexing@DNA:~/projects/zhaoxiujuan/aligned$ cat Scr_distribution.log
Total Reads                   44535835
Total Tags                    75880491
Total Assigned Tags           37717256
=====================================================================
Group               Total_bases         Tag_count           Tags/Kb
CDS_Exons           35271889            808827              22.93
5'UTR_Exons         13156148            236757              18.00
3'UTR_Exons         35031450            586598              16.74
Introns             1249039112          30136632            24.13
TSS_up_1kb          19566867            461946              23.61
TSS_up_5kb          88704854            1929741             21.75
TSS_up_10kb         163160166           3349520             20.53
TES_down_1kb        20819429            527166              25.32
TES_down_5kb        89844112            1590630             17.70
TES_down_10kb       160449429           2598922             16.20
=====================================================================

关于RSeQC的其他使用方法，参考文章：用RSeQC对比对后的转录组数据进行质控
3. 使用shell script对多个数据执行以上命令
关于脚本，参考文章：Linux学习笔记-学习shell脚本的使用（持续更新）

#脚本如下
#!/bin/bash
#program:
#  This program is running RSeQC software command and saving the output results.
#History:
#    2020/06/04            zexing         First release
for i in msh1 msh2 m3108 m3110 m3111 m3112 m3113 m3114
dobam_stat.py -i ${i}.bam.sort > ${i}_bam_stat.logread_distribution.py -i ${i}.bam.sort -r /f/xudonglab/zexing/projects/zhaoxiujuan/aligned/BED_file/hg19_RefSeq.bed > ${i}_distribution.log
done
#后台执行该命令
(base) zexing@DNA:~/projects/zhaoxiujuan/aligned$ nohup sh -x RSeQC.sh &
[1] 172776

这篇关于RNA-seq流程学习笔记（9）-使用RSeQC软件对生成的BAM文件进行质控的文章就介绍到这儿，希望我们推荐的文章对编程师们有所帮助！


http://www.chinasem.cn/article/356252。

相关文章

Python调用Orator ORM进行数据库操作

Python调用Orator ORM进行数据库操作

《Python调用OratorORM进行数据库操作》OratorORM是一个功能丰富且灵活的PythonORM库,旨在简化数据库操作,它支持多种数据库并提供了简洁且直观的API,下面我们就... 目录Orator ORM 主要特点安装使用示例总结Orator ORM 是一个功能丰富且灵活的 python O
阅读更多...
Nginx设置连接超时并进行测试的方法步骤

Nginx设置连接超时并进行测试的方法步骤

《Nginx设置连接超时并进行测试的方法步骤》在高并发场景下,如果客户端与服务器的连接长时间未响应,会占用大量的系统资源,影响其他正常请求的处理效率,为了解决这个问题,可以通过设置Nginx的连接... 目录设置连接超时目的操作步骤测试连接超时测试方法：总结：设置连接超时目的设置客户端与服务器之间的连接
阅读更多...
Java中String字符串使用避坑指南

Java中String字符串使用避坑指南

《Java中String字符串使用避坑指南》Java中的String字符串是我们日常编程中用得最多的类之一,看似简单的String使用,却隐藏着不少“坑”,如果不注意,可能会导致性能问题、意外的错误容... 目录8个避坑点如下：1. 字符串的不可变性：每次修改都创建新对象2. 使用 == 比较字符串，陷阱满
阅读更多...
Python使用国内镜像加速pip安装的方法讲解

Python使用国内镜像加速pip安装的方法讲解

《Python使用国内镜像加速pip安装的方法讲解》在Python开发中,pip是一个非常重要的工具,用于安装和管理Python的第三方库,然而,在国内使用pip安装依赖时,往往会因为网络问题而导致速... 目录一、pip 工具简介1. 什么是 pip？2. 什么是 -i 参数？二、国内镜像源的选择三、如何
阅读更多...
使用C++实现链表元素的反转

使用C++实现链表元素的反转

《使用C++实现链表元素的反转》反转链表是链表操作中一个经典的问题,也是面试中常见的考题,本文将从思路到实现一步步地讲解如何实现链表的反转,帮助初学者理解这一操作,我们将使用C++代码演示具体实现,同... 目录问题定义思路分析代码实现带头节点的链表代码讲解其他实现方式时间和空间复杂度分析总结问题定义给定
阅读更多...
Linux使用nload监控网络流量的方法

Linux使用nload监控网络流量的方法

《Linux使用nload监控网络流量的方法》Linux中的nload命令是一个用于实时监控网络流量的工具,它提供了传入和传出流量的可视化表示,帮助用户一目了然地了解网络活动,本文给大家介绍了Linu... 目录简介安装示例用法基础用法指定网络接口限制显示特定流量类型指定刷新率设置流量速率的显示单位监控多个
阅读更多...
JavaScript中的reduce方法执行过程、使用场景及进阶用法

JavaScript中的reduce方法执行过程、使用场景及进阶用法

《JavaScript中的reduce方法执行过程、使用场景及进阶用法》：本文主要介绍JavaScript中的reduce方法执行过程、使用场景及进阶用法的相关资料,reduce是JavaScri... 目录1. 什么是reduce2. reduce语法2.1 语法2.2 参数说明3. reduce执行过程
阅读更多...
如何使用Java实现请求deepseek

如何使用Java实现请求deepseek

《如何使用Java实现请求deepseek》这篇文章主要为大家详细介绍了如何使用Java实现请求deepseek功能,文中的示例代码讲解详细,感兴趣的小伙伴可以跟随小编一起学习一下... 目录1.deepseek的api创建2.Java实现请求deepseek2.1 pom文件2.2 json转化文件2.2
阅读更多...
python使用fastapi实现多语言国际化的操作指南

python使用fastapi实现多语言国际化的操作指南

《python使用fastapi实现多语言国际化的操作指南》本文介绍了使用Python和FastAPI实现多语言国际化的操作指南,包括多语言架构技术栈、翻译管理、前端本地化、语言切换机制以及常见陷阱和... 目录多语言国际化实现指南项目多语言架构技术栈目录结构翻译工作流1. 翻译数据存储2. 翻译生成脚本
阅读更多...
Linux流媒体服务器部署流程

Linux流媒体服务器部署流程

《Linux流媒体服务器部署流程》文章详细介绍了流媒体服务器的部署步骤,包括更新系统、安装依赖组件、编译安装Nginx和RTMP模块、配置Nginx和FFmpeg,以及测试流媒体服务器的搭建... 目录流媒体服务器部署部署安装1.更新系统2.安装依赖组件3.解压4.编译安装（添加RTMP和openssl模块
阅读更多...

Copyright China编程 23002807@qq.com

沪ICP备2023033072号-2