一、rRNA序列下载

在之前的教程中，我使用的NCBI的rRNA序列：https://www.jianshu.com/p/fc60dd0d0c8a

image-20220325130207993

rRNA的结构可以在KEGG Pathway数据库中查看：https://www.genome.jp/pathway/hsa03010

hsa03010

那么人究竟有多少个rRNA基因呢，在genecode页面显示：47个基因，47个转录本。

image-20220325123227719

二、运行

# 查看rRNA个数
cat human_rRNA.fasta |grep '>' |wc -l
46

# 下载下来的fa有空白行，删除fa中的空白行
sed  -i '/^$/d' human_rRNA.fasta

与genecode接近，各个数据库的统计会稍微有点出入。

构建Hisat2比对的索引

#激活小环境
conda activate rna

mkdir rRNAindex
hisat2-build -p 12 -f human_rRNA.fasta rRNAindex/human_rRNA

# 构建完成
tree
.
├── human_rRNA.fasta
└── rRNAindex    
├── human_rRNA.1.ht2    
├── human_rRNA.2.ht2    
├── human_rRNA.3.ht2    
├── human_rRNA.4.ht2    
├── human_rRNA.5.ht2    
├── human_rRNA.6.ht2    
├── human_rRNA.7.ht2    
└── human_rRNA.8.ht2

去除rRNA序列

--summary-file：输出比对结果统计文件
--no-spliced-alignment：disable spliced alignment
--no-softclip：no soft-clipping
--norc：do not align reverse-complement version of read (off)
--no-unal：不记录没比对上的reads
-p：线程数
--dta：reports alignments tailored for transcript assemblers
--un-gz：输出没有比对上的unpaired reads
-x：索引前缀

# 创建文件夹
mkdir p alignment/rRNA_dup
index_base=rRNA_1/rRNAindex/human_rRNA
outdir=alignment/rRNA_dup/
ls *gz |while read id
do
sample_name=${id%%.*}
echo "hisat2 --summary-file ${outdir}/${sample_name}_rRNA_summary.txt --no-spliced-alignment --no-softclip --norc --no-unal -p 12 --dta --un-gz ${outdir}/${sample_name}.fastq.gz -x $index_base -U ${id} | samtools view -@ 12 -Shub - | samtools sort -@ 12 -o ${outdir}/${sample_name}_rRNA_sort.bam - "
done >Filter_rRNA.sh

# 运行 qsub Filter_rRNA.sh
nohup sh Filter_rRNA.sh >Filter_rRNA.sh.log &

samtools view 参数：

-S：输入数据格式自动检测
-h：输出结果中包含表头
-u：不压缩bam文件
-b：输出 BAM文件
-@：使用线程数

Filter_rRNA.sh内容如下：

image-20220325135114401

过滤数据reads比例统计结果在*_rRNA_summary.txt文件中，样本SRR1035213_rRNA_summary.txt结果

20743537 reads; of these:
  20743537 (100.00%) were unpaired; of these:
    20742726 (100.00%) aligned 0 times
    9 (0.00%) aligned exactly 1 time
    802 (0.00%) aligned >1 times
0.00% overall alignment rate

此样本绝大部分reads都没有比对上rRNA。
小鼠数据处理同上。
后面使用过滤rRNA后的数据进行比对。