數(shù)據(jù)預(yù)處理之后是rRNA去除名惩，這一步驟是可選步驟蛮拔。
image-20220325100508016
使用Hisat2去除數(shù)據(jù)中的rRNA的原代碼，摳出來(lái)看看：
image-20220321091037228
一晃财、rRNA序列下載
在之前的教程中叨橱，我使用的NCBI的rRNA序列：http://www.reibang.com/p/fc60dd0d0c8a
image-20220325130207993
rRNA的結(jié)構(gòu)可以在KEGG Pathway數(shù)據(jù)庫(kù)中查看：https://www.genome.jp/pathway/hsa03010
hsa03010
那么人究竟有多少個(gè)rRNA基因呢，在genecode頁(yè)面顯示：47個(gè)基因断盛，47個(gè)轉(zhuǎn)錄本罗洗。
image-20220325123227719
二、運(yùn)行
# 查看rRNA個(gè)數(shù)
cat human_rRNA.fasta |grep '>' |wc -l
46
# 下載下來(lái)的fa有空白行郑临，刪除fa中的空白行
sed -i '/^$/d' human_rRNA.fasta
與genecode接近栖博，各個(gè)數(shù)據(jù)庫(kù)的統(tǒng)計(jì)會(huì)稍微有點(diǎn)出入。
構(gòu)建Hisat2比對(duì)的索引
#激活小環(huán)境
conda activate rna
mkdir rRNAindex
hisat2-build -p 12 -f human_rRNA.fasta rRNAindex/human_rRNA
# 構(gòu)建完成
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：輸出比對(duì)結(jié)果統(tǒng)計(jì)文件
- --no-spliced-alignment：disable spliced alignment
- --no-softclip：no soft-clipping
- --norc：do not align reverse-complement version of read (off)
- --no-unal：不記錄沒(méi)比對(duì)上的reads
- -p：線程數(shù)
- --dta：reports alignments tailored for transcript assemblers
- --un-gz：輸出沒(méi)有比對(duì)上的unpaired reads
- -x：索引前綴
# 創(chuàng)建文件夾
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
# 運(yùn)行 qsub Filter_rRNA.sh
nohup sh Filter_rRNA.sh >Filter_rRNA.sh.log &
samtools view 參數(shù)：
- -S：輸入數(shù)據(jù)格式自動(dòng)檢測(cè)
- -h：輸出結(jié)果中包含表頭
- -u：不壓縮bam文件
- -b：輸出 BAM文件
- -@：使用線程數(shù)
Filter_rRNA.sh內(nèi)容如下：
image-20220325135114401
過(guò)濾數(shù)據(jù)reads比例統(tǒng)計(jì)結(jié)果在*_rRNA_summary.txt文件中厢洞，樣本SRR1035213_rRNA_summary.txt結(jié)果
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都沒(méi)有比對(duì)上rRNA仇让。
小鼠數(shù)據(jù)處理同上。
后面使用過(guò)濾rRNA后的數(shù)據(jù)進(jìn)行比對(duì)躺翻。
