写在前面——之前使用的数据是单端测序,但是现在的数据基本都是双端测序。所以又找了个双端测序的例子来练习。之前在单端测序数据中,因为参考基因组注释文件找的不对,所以reads计数没有做好。这次数据质量不错,所以省略了质控和清洗,直接进入主题。由于租的服务器是2核+8G的,所以在生成sam文件和sort以及htseq-count都花费了大量的时间(四个样本集整整跑了将近一整天)。不过最后结果算是复现出来了,甚是欣慰。
文章名:AKAP95 regulates splicing through scaffolding RNAs and RNA processing factors
参考:https://www.jianshu.com/p/6d4cba26bb60
a. 建好相关文件夹
b. 00ref:存放参考基因组和基因组注释文件(红色框框内为本文需要的文件)
c. 01raw_data:双端测序,所以一个样本有两个文件。
d. clean_data:存放处理过后的数据,本文数据质量不错,所以不用清洗即可使用
e. align_data:存放比对后的文件
f. matrix:存放reads计数文件
https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE81916
具体:快速下载SRA数据
for ((i=59;i<=62;i++)) ;do prefetch -v SRR35899$i; done
fastq-dump --gzip --split-3 SRR35899*.sra
# 在另一个Linux用户下载的,需要传到自己的目录下
scp SRR35899*.gz root@dzfly:/root/project/RNA/akap95/01raw_data/
检测公司给的数据是否完整:md5sum -c md5.txt
# 使用hisat2进行比对
# -t 显示时间
# -p 设置线程
for ((i=59;i<=62;i++)); do hisat2 -t -p 2 -x 00ref/mm10/genome -1 01raw_data/SRR35899${i}_1.fastq.gz -2 01raw_data/SRR35899${i}_2.fastq.gz -S 03align_out/SRR35899${i}.sam; done
# 默认为pos排序,即默认按染色体顺序排列
for ((i=59;i<=62;i++)); do samtools sort -O bam -@ 2 -o SRR35899${i}.bam SRR35899${i}.sam; done
for ((i=59;i<=62;i++)); do samtools index SRR35899${i}.bam; done
# flagstat检查一下结果
for ((i=59;i<=62;i++)); do samtools flagstat -@ 2 SRR35899${i}.bam > SRR35899${i}.flagstat; done
具体意义例子 参考:https://cloud.tencent.com/developer/article/1703017
注意:sam文件真的很大,一个大概有30G。如果你的服务器磁盘容量不够,建议生成一个sam就将其转化为bam文件,然后将sam文件删除,再进行下一个。
for ((i=59;i<=62;i++)); do hisat2 -t -p 2 -x 00ref/mm10/genome -1 01raw_data/SRR35899${i}._1.fastq.gz -2 01raw_data/SRR35899${i}._2.fastq.gz -S 03align_out/SRR35899${i}.sam | samtools sort -O bam -@ 2 -o SRR35899${i}.bam SRR35899${i}.sam | samtools index SRR35899${i}.bam | rm SRR35899${i}.sam; done
# 方法一
# 首先将bam文件按照name进行排序,默认为pos排序
for ((i=59;i<=62;i++)); do samtools sort -@ 2 -n SRR35899${i}.bam -o SRR35899${i}_nsorted.bam; done
# 使用htseq-count进行计算
for ((i=59;i<=62;i++)); do htseq-count -r name -f bam 03align_out/SRR35899${i}_nsorted.bam 00ref/gencode.vM10.chr_patch_hapl_scaff.annotation.gtf > 04matrix/SRR35899${i}.count; done
# 方法二
# 不按name排序,直接计算
for ((i=59;i<=62;i++)); do htseq-count -r pos -f bam 03align_out/SRR35899${i}.bam 00ref/gencode.vM10.chr_patch_hapl_scaff.annotation.gtf > 04matrix/SRR35899${i}.count; done
head -n 4 SRR35899*.count
tail -n 4 SRR35899*.count
no_feature:比对区域与任何基因都没有重叠 。
ambiguous:比对区域与多个基因都发生重叠
too_low_aQual:比对质量低于设定阈值(默认是10)
not_aligned:无法比对上
alignment_not_unique:比对位置不唯一
至此,上游分析就结束了。最后获得的count文件真是来之不易呀。之后就可以使用R语言等进行分析、绘图等等了。
没有找到对应的参考基因组注释文件,导致我的htseq-count结果全都是0。后来发现参数使用的也不对。具体如下:
注意 -r name;这里我用的文件其实是按默认pos排序的bam文件,所以结果都是0!!!
htseq-count -r name -f bam 03align_out/SRR3589959.bam 00ref/Mus_musculus.GRCm38.102.gtf > 04matrix/test1.count
看一下htseq-count的具体参数 :
-f 可以自动识别文件类型
-r 默认处理按name排序的文件
这里使用 -r pos,但是结果还是0。所以是参考基因组注释文件的问题。
htseq-count -r pos -f bam 03align_out/SRR3589959.bam 00ref/Mus_musculus.GRCm38.102.gtf > 04matrix/test2.count
找到准确的参考基因组注释文件:https://www.gencodegenes.org/mouse/release_M10.html
将文件按照name排序,再进行计算,结果终于不是0了!!!
samtools sort -n SRR3589959.bam -o SRR3589959_nsorted.bam
htseq-count -r name -f bam 03align_out/SRR3589959_nsorted.bam 00ref/gencode.vM10.chr_patch_hapl_scaff.annotation.gtf > 04matrix/test3.count
这里设置 -r pos,直接对没有按name排序的文件进行计算,可以看到结果跟test3一样。( •̀ ω •́ )y
htseq-count -n 10 -r pos -f bam 03align_out/SRR3589959.bam 00ref/gencode.vM10.chr_patch_hapl_scaff.annotation.gtf > 04matrix/test4.count
这个问题告诉我们,在处理数据之前一定要确保参考基因组和注释文件是正确的。不然会浪费很多时间!!!