数据预处理之后是rRNA去除，这一步骤是可选步骤。
image-20220325100508016
使用Hisat2去除数据中的rRNA的原代码，抠出来看看：
image-20220321091037228
一、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后的数据进行比对。
