傳送門：
RNAseq005 轉(zhuǎn)錄組入門(5)：序列比對
RNAseq004 轉(zhuǎn)錄組入門(4)：參考基因組下載
RNAseq003 轉(zhuǎn)錄組入門(3)：了解fastq測序數(shù)據(jù)
RNAseq002 轉(zhuǎn)錄組入門(2)：數(shù)據(jù)下載
RNAseq001 轉(zhuǎn)錄組入門(1)：資源準(zhǔn)備
前面的五章我們分析的是人類mRNA-Seq測序的結(jié)果善涨，一般而言RNA-Seq數(shù)據(jù)分析都要有重復(fù)窒盐，而文章中有一個樣本缺少配對數(shù)據(jù)草则，所以還是選用小鼠的數(shù)據(jù)把流程再來一遍
1.數(shù)據(jù)下載及質(zhì)控見前文
2.比對
# HISAT2比對
for i in {59..62};do hisat2 -t -x /mnt/e/Work/bioinfo/public/index/mouse/hisat2/grcm38/genome -1 /mnt/e/Work/bioinfo/project/202009_RNAseq/data/SRR35899${i}_1.fastq.gz -2 /mnt/e/Work/bioinfo/project/202009_RNAseq/data/SRR35899${i}_2.fastq.gz -S /mnt/e/Work/bioinfo/project/202009_RNAseq/result/align/20200910mouse/SRR35899${i}.sam;done
3.SAM2BAM
# SAM文件轉(zhuǎn)換為BAM
for i in `seq 59 62`
do
samtools view -S SRR35899${i}.sam -b > SRR35899${i}.bam
done
4.bam flag統(tǒng)計
# 對排序后的bam統(tǒng)計flagstat
# basename命令用于獲取路徑中的文件名或路徑名，可以對末尾的擴(kuò)展名進(jìn)行刪除
ls *.bam |while read id ;do (samtools flagstat -@ 1 $id > $(basename ${id} ".bam").flagstat );done
mkdir flagstat && mv *.flagstat flagstat && cd flagstat
multiqc ./
4.1 用一個小腳本把統(tǒng)計信息轉(zhuǎn)換為csv文件
# 創(chuàng)建腳本
cat > stat.sh
### 將以下內(nèi)容寫入stat.sh
#!/bin/bash
cat *.flagstat | awk '{print $1}' | paste - - - - - - - - - - - - - > file1
# 77607517 16671207 0 0 75387881 60936310 30468155 30468155 56502696 57494864 1221810 832364 530657
# 134310379 28365145 0 0 130964009 105945234 52972617 52972617 98979648 100621038 1977826 1398380 907493
# 94264829 20737377 0 0 91921243 73527452 36763726 36763726 68525830 69723750 1460116 1023854 644490
# 111681106 24075844 0 0 109169544 87605262 43802631 43802631 82145504 83390620 1703080 1013088 643888
# 取行名
cut -d"+" -f 2 SRR3589959.flagstat | cut -d" " -f 3-90 > file2
# in total (QC-passed reads
# secondary
# supplementary
# duplicates
# mapped (97.14% : N/A)
# paired in sequencing
# read1
# read2
# properly paired (92.72% : N/A)
# with itself and mate mapped
# singletons (2.01% : N/A)
# with mate mapped to a different chr
# with mate mapped to a different chr (mapQ>=5)
# 取列名
ls *.flagstat | while read id ;do echo $(basename ${id} ".flagstat") ;done > file3
# SRR3589959
# SRR3589960
# SRR3589961
# SRR3589962
paste file3 file1 > file4
# 將file4行列轉(zhuǎn)置
awk '{
for (i=1;i<=NF;i++){
if (NR==1){
res[i]=$i
}
else{
res[i]=res[i]" "$i
}
}
}END{
for(j=1;j<=NF;j++){
print res[j]
}
}' file4 > file5
# 在file2首行加入內(nèi)容
sed '1i Index' file2 > file6
paste file6 file5 > stat.txt
cat stat.txt > stat.csv
rm file*
# 腳本內(nèi)容截止
# ==========================================================================
# 退出腳本編輯Enter登钥，ctrl+c
# 運(yùn)行腳本
bash stat.sh
4.2 csv文件處理
image.png
-
csv文件打開后是這個樣子：
image.png -
選中第一列→數(shù)據(jù)→分列→分隔符→選擇tab分隔
image.png -
選中第二列→數(shù)據(jù)→分列→分隔符→選擇空格分隔
image.png -
轉(zhuǎn)換完成
image.png
5.bam排序畔师，索引
# 排序，索引
for i in `seq 59 62`
do
samtools sort SRR35899${i}.bam -o SRR35899${i}_sorted.bam
samtools index SRR35899${i}_sorted.bam
done
# 將SAM轉(zhuǎn)換為BAM牧牢，并排序構(gòu)建索引看锉，隨后刪除SAM文件
# for i in `seq 59 62`
# do
# samtools view -S SRR35899${i}.sam -b > SRR35899${i}.bam
# samtools sort SRR35899${i}.bam -o SRR35899${i}_sorted.bam
# samtools index SRR35899${i}_sorted.bam
# done
# rm *.sam
6 注釋
# 注釋
for i in {59..62}
do
htseq-count -s no -f bam -r pos /mnt/e/Work/bioinfo/project/202009_RNAseq/result/align/20200910mouse/SRR35899${i}_sorted.bam /mnt/e/Work/bioinfo/public/Annotation/mouse/gencode/gencode.vM25.annotation.gff3 > /mnt/e/Work/bioinfo/project/202009_RNAseq/result/annotation/SRR35899${i}.count
done
# 代碼運(yùn)行報錯
# Please Install PySam to use the BAM_Reader Class (http://code.google.com/p/pysam/)Error occured when reading beginning of BAM file.
# No module named pysam
# [Exception type: ImportError, raised in __init__.py:1086]
# 解決辦法
# 下載pysam源代碼
# 下載地址：https://pypi.org/project/pysam/#files
# 復(fù)制下載鏈接放入迅雷：https://files.pythonhosted.org/packages/99/5a/fc440eb5fffb5346e61a38b49991aa552e4b8b31e8493a101d2833ed1e19/pysam-0.16.0.1.tar.gz
cd ~/biosoft
mkdir pysam && cd pysam
wget https://files.pythonhosted.org/packages/99/5a/fc440eb5fffb5346e61a38b49991aa552e4b8b31e8493a101d2833ed1e19/pysam-0.16.0.1.tar.gz
tar zxvf pysam-0.16.0.1.tar.gz
cd pysam-0.16.0.1
python setup.py install
# 報錯
# Traceback (most recent call last):
# File "setup.py", line 24, in <module>
# from setuptools import setup, find_packages
# ImportError: No module named setuptools
# python2環(huán)境下安裝setuptools
sudo apt-get install python-setuptools
# python3環(huán)境下安裝setuptools
sudo apt-get install python3-setuptools
# 再次執(zhí)行安裝
sudo python setup.py install
# 再次運(yùn)行注釋
# 構(gòu)建腳本
cat > annotation.sh
#### 輸入以下內(nèi)容
#!/bin/bash
for i in {59..62}
do
# .sorted.bam地址
input="/mnt/e/Work/bioinfo/project/202009_RNAseq/result/align/20200910mouse/SRR35899${i}_sorted.bam"
# .gtf地址
annotation="/mnt/e/Work/bioinfo/public/Annotation/mouse/gencode/gencode.vM25.annotation.gff3"
# 輸出文件地址
output="/mnt/e/Work/bioinfo/project/202009_RNAseq/result/annotation"
htseq-count -s no -f bam -r pos ${input} ${annotation} > ${output}/SRR35899${i}.count
done
# ===============================
# 運(yùn)行
bash annotation.sh
7 featureCounts統(tǒng)計
# featureCounts計數(shù)
featureCounts -p -t exon -g gene_id -a /mnt/e/Work/bioinfo/public/Annotation/mouse/gencode/gencode.vM25.annotation.gff3 -o /mnt/e/Work/bioinfo/project/202009_RNAseq/result/count/all.id.txt /mnt/e/Work/bioinfo/project/202009_RNAseq/result/align/20200910mouse/SRR35899{59..62}_sorted.bam
# 運(yùn)行后報錯
# featurecounts segmentation fault (core dumped)
# 解決辦法
# 下載二進(jìn)制版本subread
rm -rf ~/biosoft/subread
mkdir -p ~/biosoft/subread && cd ~/biosoft/subread
wget https://nchc.dl.sourceforge.net/project/subread/subread-2.0.1/subread-2.0.1-Linux-x86_64.tar.gz
tar zxvf subread-2.0.1-Linux-x86_64.tar.gz
cd subread-2.0.1-Linux-x86_64
cd ~/biosoft/subread/subread-2.0.1-Linux-x86_64/bin
./featureCounts
echo "export PATH=\$PATH:/home/cqs/biosoft/subread/subread-2.0.1-Linux-x86_64/bin" >> ~/.bashrc
source ~/.bashrc
featureCounts
# 再次運(yùn)行代碼
featureCounts -p -t exon -g gene_id -a /mnt/e/Work/bioinfo/public/Annotation/mouse/gencode/gencode.vM25.annotation.gff3 -o /mnt/e/Work/bioinfo/project/202009_RNAseq/result/count/all.id.txt /mnt/e/Work/bioinfo/project/202009_RNAseq/result/align/20200910mouse/SRR35899{59..62}_sorted.bam
# 對all.id.txt.summary進(jìn)行multiqc，查看Counts質(zhì)控
multiqc ./all.id.txt.summary
# [INFO ] multiqc : This is MultiQC v1.9
# [INFO ] multiqc : Template : default
# [INFO ] multiqc : Searching : /mnt/e/Work/bioinfo/project/202009_RNAseq/result/count/all.id.txt.summary
# Searching 1 files.. [####################################] 100%
# [INFO ] feature_counts : Found 4 reports
# [INFO ] multiqc : Compressing plot data
# [INFO ] multiqc : Report : multiqc_report.html
# [INFO ] multiqc : Data : multiqc_data
# [INFO ] multiqc : MultiQC complete
image.png
8.htseq-count統(tǒng)計
cat > htseq-count.sh
### 輸入以下內(nèi)容
#!/bin/bash
for i in {59..62}
do
# .sorted.bam地址
input="/mnt/e/Work/bioinfo/project/202009_RNAseq/result/align/20200910mouse/SRR35899${i}_sorted.bam"
# .gtf地址
annotation="/mnt/e/Work/bioinfo/public/Annotation/mouse/gencode/gencode.vM25.annotation.gff3"
# 輸出文件地址
output="/mnt/e/Work/bioinfo/project/202009_RNAseq/result/annotation"
htseq-count -s no -f bam -r pos ${input} ${annotation} > ${output}/SRR35899${i}.count
echo "SRR35899${i}.count is completed"
done
#==========================
# 運(yùn)行腳本
bash htseq-count.sh
htseq-count.sh運(yùn)行結(jié)果
