RNA-Seq学习记录（二）——入门流程命令

5软件安装在上一篇资料搜集的链接中大神们已经详细介绍

以下是我参照链接自己跑的流程命令

比对 hisat2

$ hisat2 -p 16 -x /path/index -1 R1.fq.gz -2 R2.fq.gz -S R1-2.sam

hisat2比对时报错：

File "/home/panda/miniconda3/envs/py3/bin/hisat2_read_statistics.py", line 210, in <module>

    reads_stat(args.read_file, args.read_count)

  File "/home/panda/miniconda3/envs/py3/bin/hisat2_read_statistics.py", line 168, in reads_stat

    for id, seq in fstream:

  File "/home/panda/miniconda3/envs/py3/bin/hisat2_read_statistics.py", line 44, in parser_FQ

    if line[0] == '@':

IndexError: index out of range

可参考：Question: hisat2 error index out of range ； Question: HISAT2: an "IndexError: index out of range"

samtools对SAM文件处理

sam转bam

$ samtools view -bS R1-2.SAM > R1-2.BAM

bam文件排序（-n是按照read名称，默认是-p染色体位置。这里就按照默认的染色体位置进行排序）

$ samtools sort R1-2.bam -o R1-2_sorted.bam

排序后的bam文件建立索引（按照read名称排序建立索引时报错，按照默认染色体位置时不报错，原因没有探究）

$ samtools index R1-2_sorted.bam

产生一个R1-2_sorted.bam.bai索引文件，在将R1-2_sorted.bam载入IGV可视化时可以用到！

RSeQC质控（使用方法官方网站这里，还可以参考这里和这里）

统计分析，对bam文件进行统计分析

$ bam_stat.py -i R1-2_sorted.bam

查看基因覆盖率，需要用到RefSeq.bed文件，前往RSeQC的网站这里下载，或者可以用gtf转换（还不晓得怎样转换）

$ read_distrbution.py -i R1-2_sorted.bam -r mm10_RefSeq.bed

reads计数htseq-count

这里需要下载小鼠注释文件，前往这里下载（此处GTF文件不用排序，若是传入到IGV则需要排序，可以使用IGV的tools中run igvtools的sort功能排序）

$ htseq-count -s no -r pos -f bam R1.sorted.bam gencode.vM23.chr_patch_hapl_scaff.annotation.gtf > R1.sorted_matrix.count 2> R1.sorted_matrix.log