查看原文
其他
号外:绝大部分生信技能树粉丝都没有机会加我微信,已经多次满了5000好友,所以我开通了一个微信小号,前100名添加我,仅需150元即可,3折优惠期机会不容错过哈。我的微信小号二维码在:0元,10小时教学视频直播《跟着百度李彦宏学习肿瘤基因组测序数据分析》

前些天我在生信技能树提出来了一个转录组数据分析的疑难杂症:RNA-seq的fastq文件里面为什么有gc含量的双峰,就是fastq测序数据质量控制的时候发现了GC含量的双峰,然后我简单分析了那些高重复的reads,发现一些是BCA文库,一些是rRNA相关的,所以想到了可以首先去除fastq文件中的rRNA ,就安排学徒去探索一下,过程如下

  1. 从NCBI下载rRNA的fa序列

1.1 打开NCBI,选择“Taxonomy” ,搜索 “Homo“


1.2 点击 Homo

1.3 点击 Homo sapiens


1.4 再次点击Homo sapiens


1.5 右边表格点击Nucleotide的“Subtree links”

1.6  左侧选择rRNA


1.7 下载fasta数据

注意:选择显示200个记录,让所有记录都显示在一页,方便下载全部条目。


另存为:hg38_rRNA.fasta,这个文件很小,如果大家完成上面的步骤有困难,也可以留言自己的邮箱,我们直接发送这个文件给大家哈。

  1. bowtie2建立rRNA索引

bowtie2-build hg38_rRNA.fasta hg38_rRNA

这个时候的比对工具的选择,并不一定要bowtie2软件哈。

  1. 使用bowtie2去除rRNA,重点--un-conc-gz参数。查阅hisat2的帮助文档,发现有同样的参数,所以可以用hisat2完成同样的操作。
index=/data/reference/index/bowtie2/hg38_rRNA/hg38_rRNA
bowtie2 -x ${index} --un-conc-gz SRR6236728_rmrRNA.fq.gz -1 SRR6236728_1_val_1.fq.gz -2 SRR6236728_2_val_2.fq.gz -p 6 | samtools sort -o tmp.bam -
rm tmp.bam

bowtie2运行耗时20min左右,如果是hisat2,可能会更快,大家可以自行测试和比较一下。

  1. 去除前后的fq文件比较


上面fq文件包含val标签的是去除rRNA之前的fq文件,包含rmrRNA标签的是去除rRNA之后的fq文件。


下面是不同fq文件的fastqc报告:

可以看到,少了8百万条reads。



可以看到,GC含量最后一个峰是由rRNA导致,因为8百万条reads被去除后,该峰就消失了。

最后剩下的问题,就是GC含量的另外一个峰。我们后续再谈它!


由于最近订阅号消息中,微信文章不再按照时间顺序排列,导致你可能无法及时看到我们的生信技能树的教程,所以我想邀请你:将生信技能树设为“星标”或经常为文章点“在看”。就是这里👇

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

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