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

2023-11-06 11:10

本文主要是介绍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构建智能BAT文件生成器的完美解决方案

《使用Python构建智能BAT文件生成器的完美解决方案》这篇文章主要为大家详细介绍了如何使用wxPython构建一个智能的BAT文件生成器,它不仅能够为Python脚本生成启动脚本,还提供了完整的文... 目录引言运行效果图项目背景与需求分析核心需求技术选型核心功能实现1. 数据库设计2. 界面布局设计3

使用IDEA部署Docker应用指南分享

《使用IDEA部署Docker应用指南分享》本文介绍了使用IDEA部署Docker应用的四步流程:创建Dockerfile、配置IDEADocker连接、设置运行调试环境、构建运行镜像,并强调需准备本... 目录一、创建 dockerfile 配置文件二、配置 IDEA 的 Docker 连接三、配置 Do

Android Paging 分页加载库使用实践

《AndroidPaging分页加载库使用实践》AndroidPaging库是Jetpack组件的一部分,它提供了一套完整的解决方案来处理大型数据集的分页加载,本文将深入探讨Paging库... 目录前言一、Paging 库概述二、Paging 3 核心组件1. PagingSource2. Pager3.

Python进行JSON和Excel文件转换处理指南

《Python进行JSON和Excel文件转换处理指南》在数据交换与系统集成中,JSON与Excel是两种极为常见的数据格式,本文将介绍如何使用Python实现将JSON转换为格式化的Excel文件,... 目录将 jsON 导入为格式化 Excel将 Excel 导出为结构化 JSON处理嵌套 JSON:

python使用try函数详解

《python使用try函数详解》Pythontry语句用于异常处理,支持捕获特定/多种异常、else/final子句确保资源释放,结合with语句自动清理,可自定义异常及嵌套结构,灵活应对错误场景... 目录try 函数的基本语法捕获特定异常捕获多个异常使用 else 子句使用 finally 子句捕获所

C++11右值引用与Lambda表达式的使用

《C++11右值引用与Lambda表达式的使用》C++11引入右值引用,实现移动语义提升性能,支持资源转移与完美转发;同时引入Lambda表达式,简化匿名函数定义,通过捕获列表和参数列表灵活处理变量... 目录C++11新特性右值引用和移动语义左值 / 右值常见的左值和右值移动语义移动构造函数移动复制运算符

Python对接支付宝支付之使用AliPay实现的详细操作指南

《Python对接支付宝支付之使用AliPay实现的详细操作指南》支付宝没有提供PythonSDK,但是强大的github就有提供python-alipay-sdk,封装里很多复杂操作,使用这个我们就... 目录一、引言二、准备工作2.1 支付宝开放平台入驻与应用创建2.2 密钥生成与配置2.3 安装ali

C#中lock关键字的使用小结

《C#中lock关键字的使用小结》在C#中,lock关键字用于确保当一个线程位于给定实例的代码块中时,其他线程无法访问同一实例的该代码块,下面就来介绍一下lock关键字的使用... 目录使用方式工作原理注意事项示例代码为什么不能lock值类型在C#中,lock关键字用于确保当一个线程位于给定实例的代码块中时

MySQL 强制使用特定索引的操作

《MySQL强制使用特定索引的操作》MySQL可通过FORCEINDEX、USEINDEX等语法强制查询使用特定索引,但优化器可能不采纳,需结合EXPLAIN分析执行计划,避免性能下降,注意版本差异... 目录1. 使用FORCE INDEX语法2. 使用USE INDEX语法3. 使用IGNORE IND

C# $字符串插值的使用

《C#$字符串插值的使用》本文介绍了C#中的字符串插值功能,详细介绍了使用$符号的实现方式,文中通过示例代码介绍的非常详细,需要的朋友们下面随着小编来一起学习学习吧... 目录$ 字符使用方式创建内插字符串包含不同的数据类型控制内插表达式的格式控制内插表达式的对齐方式内插表达式中使用转义序列内插表达式中使用