转录组RNA-SEQ上游分析2020

安装配置conda

使用清华源下载sh脚本并安装

# 使用清华源下载sh脚本
wget -c  https://mirrors.tuna.tsinghua.edu.cn/anaconda/miniconda/Miniconda3-latest-Linux-x86_64.sh

# 从官网下载最新版Miniconda3安装包，但速度较慢
wget -c https://repo.anaconda.com/miniconda/Miniconda3-latest-Linux-x86_64.sh

下载完成后直接运行脚本文件bash Miniconda3-latest-Linux-x86_64.sh。需要输入yes然后等待安装完毕
最后安装好后，还不能马上使用conda，需要source一下bashrc

# 激活bashrc
source ~/.bashrc

==注意⚠️：==

  • conda会在bashrc中写入脚本，连接ssh自动进入conda环境的命令。如果不需要可以运行命令及性能配置conda config --set auto_activate_base false
  • 另外如果使用zsh等工具如果没有自动写入zshrc，可以在文件中手动写入。
  • 如果conda命令不被读取，可以手动定义环境变量export PATH="/home/super/miniconda3/bin:$PATH"

设置镜像源

# 下面这四行配置清华大学的bioconda的channel地址，国内用户推荐
conda config --add channels https://mirrors.tuna.tsinghua.edu.cn/anaconda/pkgs/free/
conda config --add channels https://mirrors.tuna.tsinghua.edu.cn/anaconda/pkgs/main/
conda config --add channels https://mirrors.tuna.tsinghua.edu.cn/anaconda/cloud/conda-forge/
conda config --add channels https://mirrors.tuna.tsinghua.edu.cn/anaconda/cloud/bioconda/

## 官网默认
conda config --add channels r 
conda config --add channels conda-forge 
conda config --add channels bioconda

在设置后镜像或者设置不自动进入base后，会在.condarc文件中自动生成config信息。如下：

$ cat .condarc 

auto_activate_base: false
channels:
  - https://mirrors.tuna.tsinghua.edu.cn/anaconda/cloud/bioconda/
  - https://mirrors.tuna.tsinghua.edu.cn/anaconda/cloud/conda-forge/
  - https://mirrors.tuna.tsinghua.edu.cn/anaconda/pkgs/main/
  - https://mirrors.tuna.tsinghua.edu.cn/anaconda/pkgs/free/
  - defaults

conda环境创建

创建一个python2的环境管理：

conda create -y -n rna_seq python=3

# -y        自动确认
# -n        新环境名字
# python=3  新环境中python=3

激活和退出环境

conda activate <conda_name>     #激活某环境
conda decativate <conda>        #取消激活某环境

conda安装软件

在软件环境中使用命令安装软件

conda install -y sra-tools      #安装sra-tool软件，可以通过空格安装多个软件
conda install -y sra-tools fastqc trim-galore hisat2 subread multiqc samtools salmon fastp

conda软件安装位置和普通软件安装位置不一样，通过which <softname>来查看conda安装的软件位置

质量评估 @ ==fastQC==

fastq格式

FastQ格式描述：https://mp.weixin.qq.com/s/8g-oUjiEhV4cGMJNuhmISQ
FastQ格式wiki：https://en.wikipedia.org/wiki/FASTQ_format
FastQ格式文献：https://www.ncbi.nlm.nih.gov/pmc/articles/PMC2847217/

概念
FastQ格式是序列格式中常见的一种，它存储了生物序列以及相应的质量评价，其序列以及质量信息都是使用一个ASCII字符标示，最初由Sanger开发，目的是将FASTA序列与质量数据放到一起，目前已经成为高通量测序结果的事实标准。

格式说明
FASTQ文件中每个序列通常有四行：

  • 1.第一行：必须以“@”开头，后面跟着唯一的序列ID标识符，然后跟着可选的序列描述内容，标识符与描述内容用空格分开；
  • 2.第二行：序列字符（核酸为[AGCTN]+，蛋白为氨基酸字符）；
  • 3.第三行：必须以“+”开头，后面跟着可选的ID标识符和可选的描述内容，如果“+”后面有内容，该内容必须与第一行“@”后的内容相同；
  • 4.第四行：碱基质量字符，每个字符对应第二行相应位置碱基或氨基酸的质量，该字符可以按一定规则转换为碱基质量得分，碱基质量得分可以反映该碱基的错误率。这一行的字符数与第二行中的字符数必须相同。

FsatQC软件

FastQC质量评估软件官网：http://www.bioinformatics.babraham.ac.uk/projects/fastqc/
==注意⚠️==

  • fastqc可以对 *.bam *.sam *.fq *.fq.gz进行质量评估。
  • fastqc可以通过-t指定多线程操作，多线程是同时处理多个输入文件，几个线程可以同时处理几个文件，单个文件使用多线程似乎没有意义
  • 对bam质量评估和对过滤后、指控后的文件使用fastqc似乎没有区别
  • 在bash中批处理比较简单，但是zsh中，不太一样，需要在命令替换出使用 echo $list

常用参数：

# 常用参数
fastqc -o <out.dir> -t <thred_num> -f <input_format>  <input_file_1> <input_file_2> ...

# -o    设置输出目录
# -t    设置线程数
# -f    设置输入文件格式

批处理

# bash中
a=`ls *.fq`
fastqc -o ./fastqc_raw -t 10 $a

# zsh中
b=`ls -C *.fq`
fastqc -o ./fastqc_raw -t 10 `echo $b`

multiqc综合所有qc

使用multiqc来把所有的质量评估放在一起观察

multqc -o <out.path> *.fastqc.zip

结果解读

单一碱基占比

image

unique reads 的占比，

单碱基测序质量在150bp上的分布情况

image

也有这样的

image

当前RNA-seq测序技术，测序错误率分布存在以下两个特征。

  1. 测序错误率随着测序序列(Sequenced Reads)长度的增加而升高。这是由测序过程中化学试剂的消耗导致的，为Illumina高通量测序平台所具有的特征。
  2. 前6个碱基具有较高的测序错误率，此长度恰好为RNA-seq建库过程中反转录所需的随机引物长度。前6个碱基测序错误率较高是因为随机引物和RNA模版的不完全结合(Jiang et al.)。

测序质量的分布图

image

和单碱基测序质量在150bp上的分布情况不同，这个是单个样本中，碱基质量的分布情况。绝大部份集中在Q30以上，效果良好（类似于密度分布图）

GC含量测定

1.整体GC含量测定，主要看是否有双峰，如果有双峰可能有其他物种掺入。（GC含量在物种间存在一定特异性）


image

2.单碱基GC(TAN)在150bp上的分布情况。理想情况是四条线在25%轻微波动，但是如所见，前几个bp非常不稳定。这是由于反转录过程中所使用的6bp随机引物，会引起前几位碱基在核苷酸组成上有一定偏好性，产生正常波动，随后则趋于稳定。

image

image

下图除了加入了ATCG等碱基，还加入了不确定碱基N碱基的含量。后期会去除这些含N碱基的reads。

接头adapter统计

image

统计了接头出现的累积频率。一般累积频率不超过5%即可

过滤reads @ ==trim_galore==

过滤标准

一般的过滤标准：

  • 去除含有接头的reads
  • 过滤去除低质量值的reads
  • 去除含有N比例大于5%的reas

trim_galore

trim_galore使用手册：https://github.com/FelixKrueger/TrimGalore/blob/master/Docs/Trim_Galore_User_Guide.md
阅读延伸：https://www.jianshu.com/p/7a3de6b8e503

Trim Galore是对FastQC和Cutadapt的包装。适用于所有高通量测序，包括RRBS(Reduced Representation Bisulfite-Seq ), Illumina、Nextera 和smallRNA测序平台的双端和单端数据。主要功能包括两步：
第一步首先去除低质量碱基，然后去除3' 末端的adapter, 如果没有指定具体的adapter，程序会自动检测前1million的序列，然后对比前12-13bp的序列是否符合以下类型的adapter:

==注意⚠️==

  • 在conda中，trim_galore叫做trim-galore
  • length选项不要设置太高
  • -j\--cores不要设置超过4，因为设置为4实际使用内核是15
  • trim_galore似乎需要一个==python3==的环境。因此所建立的conda环境应该是python3，如果当前环境是python2，那就新建一个python3环境来运行
  • fq.gz可以不用解压也可以运行
# 常用参数
trim_galore

-j                  #使用线程数,注意假设已使用Python 3并安装了Pigz，那么内核设置为4，实际使用内核是15，因此，最高设为4
-q                  #切除质量得分低于设定值的序列
--phred33           #得分标准（默认）
-a                  #输入adapter序列，也可以不输入，软件自动搜索可能性最高的平台对应的adapter，也可以直接设置参数输入这三个平台分别是--illumina --nextera --small_ma
--illumina          #设置测序平台是ilumina，用来设定接头序列用
--length            #设置小于此bp的reads会被删除，注意，在pe150下，不要设置太高，可以50或36（默认20）
--max_length        # 设置长度大于此值被丢弃
-s/--stringency     #设定可以容忍的前后adapter的重叠碱基数，默认是1，可以放宽以为后一个adapter几乎测不到
--paired            #上端reads，如果一个被剔除，则另一个也被剔除
-o                  #设定已经存在的输出目录
--fastqc            #当分析结束后，使用默认选项对结果文件进行fastqc分析
--gzip              #输出结果使用gzip压缩,gzip压缩将大大减慢多核进程，因此几乎不值得
--max_n             #设置如果超过设定n就删除。可以根据length选项，
--basename          #设置输出文件的基本名称，尽在使用1个文件或2个文件参数时候可用

例子

trim_galore --phred33 -q 20 --length 36 --stringency 3 --fastqc --paired --max_n 3 -j 4 -o ./trim_galore 1.fastq.gz 2.fastq.gz

批处理

ls *1.fq >1.txt
ls *2.fq >2.txt
paste 1.txt 2.txt >config

# cat trim_galore.sh
cat config |while read id
do
    arr=(${id})
    fq1=${arr[0]}
    fq2=${arr[1]}
    echo "trim_galore work_ID：${fq1} started at $(date)"
    echo "trim_galore -q 20 --phred33 --stringency 3 --length 36 --paired --fastqc  --max_n 3 -j 4 ${fq1} ${fq2}"
    echo "trim_galore work_ID：${fq1} finished at $(date)"
    echo -e "\n"
done

# 进行上述脚本
nohup sh trim_galore.sh >trim_galore.log &

批处理参考：https://www.jianshu.com/p/5a77b4156a5d
并行处理要求：https://www.jianshu.com/p/74a328f4d23a

比对 @ ==hisat2==

hisat用户手册：http://ccb.jhu.edu/software/hisat2/manual.shtml
==简书：有关hisat深入介绍，推荐==：https://www.jianshu.com/p/ce3f4afb9b60

hisat的优势：

  • 更适合RNA-seq，因为索引策略的使用，可以相对轻松比对跨区域的read（可变剪切）
  • Tophat2耗时久，STAR吃内存。hisat兼有两个的耗时、内存消耗的优点

参考基因组，gtf文件选择，不同基因组版本差异

==注意⚠️：==

  • 参考基因组选择正确的版本还是有必要的。primary的版本中是不包括haplotype info的，而top level中会包含大量的变异信息，而这部分是很冗余并且一般也用不太到。primary_assembly 和 toplevel相比不包含haplotype, 更适合用于比对。
  • 解压参考基因组，会非常大
  • 对于mask/unmask 通常选择softmask或者unmasked, 一般不用rm的
  • 小鼠的基因是Mus musculus，例如ENSEMBL

参考基因组和gtf文件可以从esemble下载。资源汇总地址：ftp://ftp.ensembl.org/pub/ 。可以选择最新的 101版本（截止2020/08/25）。

image

在101版本中选gtf来下载最新的gtf注释文件，选择fasta来选择最新的基因组文件。
基因组文件依次选择release-101 > fasta > homo_sapiens > dna 然后在诸多文件中选择*GRCH38.dna.primary_assembly*
image

不同基因组版本差异：

  • 一般是选择primary来比对
  • 不同mask的差别：
    • dna 版本对应测序组装得到的原始基因组序列
    • dna_sm版本对应 用小写字母代替了重复和低复杂度区域的全基因组序列
    • dna_rm版本对应 用 “NNN” 屏蔽了重复和低复杂度区域的全基因组序列
    • 三种全基因组在序列长度上是一致的，只是在重复序列区或低复杂度区存在差异。事实上很多比对软件对大小写不敏感，例如BWA，所以用 unmasked 和 softmasked 的基因组做参考序列，比对得到的结果有时候是一致的。

参考&资料：

构建索引

人的基因组索引可以下，但是可能会有基因组差异问题，所以，除非非常清楚和索引的基因组匹配，否则建议自己构建索引。当然，构建索引是一个非常漫长的过程。

hisat的索引策略

  • hisat该软件应用了两类不同的索引类型：代表全基因组的全局FM索引和大量的局部小索引，每个索引代表64000bp
  • 针对人类基因组创建了48000个局部索引，每一个覆盖1024bp，最终可以覆盖这个3 billion 的碱基的基因组。这种存在交叉（overlap）的边界可以轻松的比对那些跨区域的read（可变剪切体）

hisat的比对策略

  • RNA-seq的read大概分为五类：
    • 不跨外显子的read
    • 长锚定read，两个外显子上都至少有16bp，约占25%，大多数可以被唯一定位到基因组
    • 中锚定read，两个外显子中，最短的有8-15bp，约占5%，利用局部索引可以实现快速检索
    • 短锚定read，两个外显子中，最短的有1-7bp，约4.2%，因为序列短，因此只能通过在hisat比对其它read时发现的剪切位点或者用户自己提供的剪切位点来辅助比对
    • 跨多外显子read，跨三个或三个以上read
image

索引的构建

注意，如果要用到gtf中exon和split site信息，要先用hisat2自带的py软件生成，然后添加到参数--exon 和参数--ss后

hisat2-build [options] <reference_in> <ht2_index_base>
# 常用[options]：
# -p                线程数
# 参数：
# <reference_in>    参考基因组fasta文件 list，如果为list，使用逗号分开
# <ht2_base>        索引文件的前缀名
# --ss              > 与 --exon 一起使用，提供拼接位点信息。只支持hisat2自己格式
#                   ^（需要用hisat2_extract_splice_sites.py 与 gtf来生成）
# --exon            指定外显子列表（只支持hisat2自己格式，需要通过hisat2_extract_exons.py 与gtf文件来生成）从GTF提取剪接位点。


# 例子
nohup \ 
hisat2-build -p 6 Homo_sapiens.GRCh38.dna.primary_assembly.fa homo_GRCh38.dna.primary \
&

索引构建完毕如下：

image

hisat2比对

# 常用参数
hisat2 [option] -1 <m1> -2 <m2> -x <index_base> -S <out.sam>

-x <path/to/base>             #索引的前缀,也可以添加路径，例如 ～/genome/homo_GRCh38.dna.primary
-1              #双端测序的第一个文件，若有多组，则用逗号隔开，名字与2要匹配，如-1 flyA_1.fq,flyB_1.fq
-2              #类比上，名字与1要匹配，如-2 flyA_2.fq,flyB_2.fq
-U              #单端数据文件
-S              #保存到的sam文件，默认输出到标准输出

# 以下为参数
-p              #线程数
--dta           # > 在下游使用stringtie组装的时候量身定做的参数。使用此选项，
                # ^ HISAT2需要更长的锚长度才能从头发现拼接位点，这样可以减少与短锚的对齐，
                # ^ 这有助于转录汇编程序显着提高计算和内存使用率。
                # ^ 当然下游不使用stringtie也可以使用此参数减少短锚定read
--dta-cufflinks #下游使用cufflinks则需要添加此参数
-q              #指定读取的测序文件是fastq文件（含有质量信息）
-f              #指定读取的测序文件是fa文件
-t              #打印加载索引文件和对齐读数所需的时钟时间
--phred33       #质量表示方法，如果使用fq文件，建议选择，同理有小众的--phred64
--min-intronlen #设置最小内含子的长度，默认值20
--max-intronlen #设置最大内含子长度，默认500000
--known-splicesite-infile <path>    #提供已知的剪接位点列表，HISAT2利用这些列表比对,可以用hisat2_extract_splice_sites.py和gtf生成信息
--novel-splicesite-outfile <path>   #在结果中生成一个剪接位点的列表
--novel-splicesite-infile <path>    # > 可以使用novel-splicesite-outfile所生成的剪接列表作为
                                    # ^ 新剪接列表，应该可以自己定义。为特定基因。

例子

# 例子 
hisat2 -p 10 -q --known-splicesite-infile <splicesites.txt> -x <homo_GRCh38_dna_primary> -1 <S_CTRL_1_val_1.fq.gz>,<S_OE_1_val_1.fq.gz> -2 <S_CTRL_2_val_2.fq.gz>,<S_OE_2_val_2.fq.gz> -S <./sam/sample.hisat.sam>

批处理

ls *1.fq.gz|cut -d"_" -f1,2>id.txt      #整理列举id，在trim_galore目录下，用已经过滤好的数据

# cat hista_align.sh
cat ~/lab/ZhangJinRui/analysis/trim_galore/id.txt | while read id                                
do
    echo "hista2 :work_ID ${id} start at $(date)"
    echo "hisat2 -p 10 -q --known-splicesite-infile ~/lab/ZhangJinRui/analysis/geom/splicesites.txt -x homo_GRCh38_dna_primary -1 ~/lab/ZhangJinRui/analysis/trim_galore/${id}_1_val_1.fq.gz -2 ~/lab/ZhangJinRui/analysis/trim_galore/${id}_2_val_2.fq.gz -S ~/lab/ZhangJinRui/analysis/sam/${id}.hisat.sam"
    echo "hisat2 :work_ID ${id} finished at $(date)"
done

# 运行上述脚本 在索引目录下？不用也可以，索引前可以用路径指定 -x 
nohup sh hista_align.sh >hisat_align.log &

PS：
使用hisat2_extract_splice_sites.py来生成gtf中的剪接位点，辅助比对。使用方法：

hisat2_extract_splice_sites.py genes.gtf > splicesites.txt

sam格式转换 @ ==samtools==

sam格式：The Sequence Alignment / Map format，存储了序列比对情况的文件，是一个文本文件，都以tab分列。体积较大。分为头部区和主体区。
官方文档：http://samtools.github.io/hts-specs/SAMv1.pdf

头部区

image
  • @HD VN:1.0 SO:unsorted：sam文件的第一行
    • VN 表示是格式版本
    • SO 表示排序类型，有unknown（default），unsorted，queryname和coordinate（一般要安坐标排序）
  • @SQ SN:1 LN:248956422
    • SN 表示参考序列名，这里面是1号染色体，一般是1-22 + X+Y 染色体，也可能由于使用基因组的不同导致有其他
    • LN 是该参考序列长度，此处指1号染色体长度
  • @PG ID:hisat2 PN:hisat2 VN:2.2.1 CL:"/home/super/miniconda3/envs/rna_seq3/bin/hisat2-align-s --wrapper basic-0 -p 10 .......
    • 软件相关信息

其他具体sam格式信息：https://blog.csdn.net/genome_denovo/article/details/78712972

PS:诺和分析数据的代码：

/TJPROJ6/RNA_T/software/miniconda2/envs/hisat2-2.0.5/bin/hisat2-align-s --wrapper basic-0 -x /TJPROJ6/GB_TR/reference_data/new_pip/Animal/Homo_sapiens/Homo_sapiens_Ensemble_94/Homo_sapiens_Ensemble_94 -p 4 --dta -t --phred33 --passthrough -1 /tmp/22775.inpipe1 -2 /tmp/22775.inpipe2

排序和转BAM

samtools官方文档：http://www.htslib.org/doc/samtools.html
==注意⚠️：==

  • samtools默认使用坐标排序
  • 当未指定samtool的索引类型时，默认是使用BAI索引。也可以samtools -c 来指定CSI索引

排序转bam

详细使用方法可以使用man samtools sort查看

samtools sort [option] <input.sam>
            -@                      #线程数
            -o <path/to/out.bam>    #输出的sort后文件
            -T <path/to/temp>       #临时文件
            -O <sam|bam..>          #输出文件格式，默认根据-o来确定
            -n                      #按照QNAME字段排序而不是染色体坐标
            
# 例子
samtools sort -@ 10 -o ./sort_bam/sample_sort.bam sample.sam

批处理

ls *.sam|cut -d"." -f1,2 >sam.id.txt        #把诸如 S_CTRL.hisat.sam S_OE.hisat.sam 取前缀名然后分行存储在文本内

# cat sam_sort.sh 
cat sam.id.txt|while read id  
do
    echo -e "*******WORK_ID ${id} START AT : \n$(date)"
    samtools sort -@ 10 -o ./sort_bam/${id}.sort.bam ${id}.sam
    echo -e "*******WORK_ID ${id} FINISHED AT : \n$(date)"
done

# 运行脚本 在当前sam文件目录下，并且在目录下新建了 sort_bam文件
nohup sh sam_sort.sh >sam_sort.log &

索引bam文件

samtools index sample.sort.bam  #生成文件默认是原文件名 + .bai 也可以指定在原文件后面直接指定

批处理

ls *.bam > bam.list.txt

#cat bam_index.sh
cat ba.list.txt|while read id
do
    echo -e "*******WORK_ID ${id} START AT : \n$(date)"
    samtools index ${id}
    echo -e "*******WORK_ID ${id} FINISHED AT : \n$(date)"
done 

#运行
nohup sh bam_index.sh >bam_index.log &

==注意⚠️：==
说明网站上指出，samtools index 有一个 -@ 指定线程的参数。但是自己实际操作报错。

计数 @ ==featureCounts==

英文简明说明书：http://bioinf.wehi.edu.au/featureCounts/
官方网站：http://subread.sourceforge.net/
英文用户手册下载：http://bioinf.wehi.edu.au/subread-package/SubreadUsersGuide.pdf （featurecounts说明在P31）
简书详细介绍（推荐）：https://www.jianshu.com/p/9cc4e8657d62

featureCounts理解的核心是feature。我们感兴趣的基因组特征（feature）可能是外显子、基因、启动子区域、gene body或者是基因组间隙，不单单是gene_id。

在计数中RNA-seq计数可能更加复杂，：

  • 因为需要考虑到外显子剪切。一种计数方法是数一下与每一个被注释的外显子重合的read，另外一种方法是数一下与每一个基因区域重合的read
  • 另一个问题是，尽管read count中包含感兴趣基因区域的所有信息，但是我们无法区分isoform的信息，因为同一基因下的不同的isoform存在很大程度的重合区域，很多基于模型的方法被开发同时也有人统计过isoform中不重合的区域作为统计的依据。

输入的GTF的作用：
GFF/GTF/SAF主要提供feature identifier(如geneID), chromosomename, start position, end position and strand 这五列信息

gtf信息解读：
https://www.jianshu.com/p/7cec5a02768a

feature和meta-feature：

  • feature是指基因组上被定义的一个片段区域，meta-feature是指多个feature组成的区域，如exon和gene的关系
  • featurecounts可以对features和meta-feature进行计数

read和fragment：
如果是双端测序，这一对read定义了一个DNA/RNA片段的两端，这种情况下，featurecounts会计算片段数（fragment）而不是read数。即一个fragment由两个raed确定。（？fragment可能包含有read中没有的bp信息，但是read负责确定fragment的定位信息）

多重overlap的选择：

  • 多重overlap是指一个read/fragment跨越了两个feature或在计数meta-feature时跨越两个meta-feature
  • 可以自己定义怎么处理这种read/fragment，如果是RNA-seq实验，我们推荐把这种read去掉

==注意⚠️：==

  • featureCounts已经整合到subread软件中，使用featureCounts可以下载subread然后就可以直接使用了
  • SAM/BAM文件和GTF/GFF/SAF文件需要来自同一个参考基因组，即必须参考基因组和GTF/GFF/SAF文件来自同一个网站，同一个版本
  • 如果是RNA-seq，则最好不要允许read/fragment的overlap(不要加 -O 参数)，如果是chip-seq实验，我们建议保留
  • 如果在计数meta-feature的时候，在同一个meta-feature中overlap两个feature的read/fragment只计数一次
  • 当想对feature水平进行统计的时候，需要设置-f参数，
  • 想对meta-feature水平进行统计的时候，不能设置-f参数
  • -g参数用来设定meta-feature；默认为gene_id，可以选择transcript_id、gene_name、transcript_name等，具体可以去GTF文件中查看属性列
  • 对于count结果，如果想了解每个基因上的count数，则只需要提取出第1列和第7列的信息
  • 数据深度不够可以不用-B剔除数据
  • 实际数据，在RNA-seq结果中，加上-B -C的筛选选项，会比不加少约5% reads。

主要参数

featureCounts [options] <input.file>

     <input.file>       #输入文件，sam/bam 支持多个文件输入
     -a <gtf>           #参考gtf注释文件，支持gz压缩
     -F <参考文件格式>    #默认GTF
     -T <int>           #线程数 1-32
     -J                 #对可变剪接计数
     -G <str>           #有-J设置时，用来提供参考基因组辅助寻找可变剪接
     -M                 #如果设置了-M则多重比对会被统计
     -o <out.file>      #输出文件的名字，输出文件的内容是read的统计数目
     -O                 #允许多重比对，当一个read比对到多个feature或多个metafeature的时候，这条read会被统计多次
     -p                 #声明测序类型是paired-end，此时，后续统计fragment而不是read
     -B                 #在-p设置时，只有两端都比对上的fragment才会被统计
     -C                 #在-p设置时，-C声明比对到不同染色体的fragment不计数
     -d <int>           #最短的fragment ，默认50
     -D <int>           #最长的fragment，默认600
     -f                 #设置后会统计feature水平数据，如exon-level，否则会统计meta-feature层面数据，如gene-level，（？应该是和-g冲突,一般与-t连用?）
     -g <str>           #参考的gtf提供的时候，我们需要提供一个id identifier 来将feature水平的统计汇总为meta-feature水平的统计，默认为gene_id，可以选择transcript_id、gene_name、transcript_name等
     -t <str>           #设置feature-type，-t指定的必须是gtf中有的feature，同时read只有落到这些feature上才会被统计到，默认是“exon”（？应该和-g冲突）

例子：

  • 普通例子
featureCounts -T 10 -a $gtf -o all.id.txt -p -t exon -g gene_id   *.sorted.bam
  • 多个bam文件
featureCounts -t exon -g gene_id -a annotation.gtf -o counts.txt library1.bam library2.bam library3.bam
  • 双端测序数据
featureCounts -p -t exon -g gene_id -a annotation.gtf -o counts.txt mapping_results_PE.bam
  • 排除映射到不同染色体的fragment/read
featureCounts -p -C -t exon -g gene_id -a annotation.gtf -o counts.txt mapping_results_PE.bam
  • 只计数那些双端都匹配的reads
featureCounts -p -B -t exon -g gene_id -a annotation.gtf -o counts.txt mapping_results_PE.bam
  • 个人例子
#输出gene_id
nohup featureCounts -T 10 -a ~/lab/ZhangJinRui/analysis/geom/Homo_sapiens.GRCh38.101.gtf -o counts.txt -p -t exon -g gene_id ~/lab/ZhangJinRui/analysis/sam/sort_bam/*.sort.bam &

#输出gene_name
nohup featureCounts -T 10 -a ~/lab/ZhangJinRui/analysis/geom/Homo_sapiens.GRCh38.101.gtf -o gene_name_counts.txt -p -t exon -C -B -g gene_name ~/lab/ZhangJinRui/analysis/sam/sort_bam/*.sort.bam >featurecount_gene_name.log &

结果处理

一般输出结果大于7列，其中第一列是由-g决定的meta-feature，第七列开始是对应bam文件的count数，第七列后面多样本的couts数，如果要引用多个数，例如，共20个样本，要选从第七列选到26列，可以使用

cat counts.txt|head|cut -f1,`seq -s"," 7 26`

我想建立并管理一个高质量的生信&统计相关的微信讨论群，如果你想参与讨论，可以添加微信：veryqun 。我会拉你进群，当然有问题也可以微信咨询我。










最后编辑于 ：
©著作权归作者所有,转载或内容合作请联系作者
  • 人面猴
    序言：七十年代末，一起剥皮案震惊了整个滨河市，随后出现的几起案子，更是在滨河造成了极大的恐慌，老刑警刘岩，带你破解...
    沈念sama阅读 194,088评论 5赞 459
  • 死咒
    序言：滨河连续发生了三起死亡事件，死亡现场离奇诡异，居然都是意外死亡，警方通过查阅死者的电脑和手机，发现死者居然都...
    沈念sama阅读 81,715评论 2赞 371
  • 救了他两次的神仙让他今天三更去死
    文/潘晓璐 我一进店门，熙熙楼的掌柜王于贵愁眉苦脸地迎上来，“玉大人，你说我怎么就摊上这事。” “怎么了？”我有些...
    开封第一讲书人阅读 141,361评论 0赞 319
  • 道士缉凶录：失踪的卖姜人
    文/不坏的土叔 我叫张陵，是天一观的道长。 经常有香客问我，道长，这世上最难降的妖魔是什么？ 我笑而不...
    开封第一讲书人阅读 52,099评论 1赞 263
  • 港岛之恋（遗憾婚礼）
    正文 为了忘掉前任，我火速办了婚礼，结果婚礼上，老公的妹妹穿的比我还像新娘。我一直安慰自己，他们只是感情好，可当我...
    茶点故事阅读 60,987评论 4赞 355
  • 恶毒庶女顶嫁案：这布局不是一般人想出来的
    文/花漫 我一把揭开白布。 她就那样静静地躺着，像睡着了一般。 火红的嫁衣衬着肌肤如雪。 梳的纹丝不乱的头发上，一...
    开封第一讲书人阅读 46,063评论 1赞 272
  • 城市分裂传说
    那天，我揣着相机与录音，去河边找鬼。 笑死，一个胖子当着我的面吹牛，可吹牛的内容都是我干的。 我是一名探鬼主播，决...
    沈念sama阅读 36,486评论 3赞 381
  • 双鸳鸯连环套：你想象不到人心有多黑
    文/苍兰香墨 我猛地睁开眼，长吁一口气：“原来是场噩梦啊……” “哼！你这毒妇竟也来了？” 一声冷哼从身侧响起，我...
    开封第一讲书人阅读 35,175评论 0赞 253
  • 万荣杀人案实录
    序言：老挝万荣一对情侣失踪，失踪者是张志新（化名）和其女友刘颖，没想到半个月后，有当地人在树林里发现了一具尸体，经...
    沈念sama阅读 39,440评论 1赞 290
  • 护林员之死
    正文 独居荒郊野岭守林人离奇死亡，尸身上长有42处带血的脓包…… 初始之章·张勋 以下内容为张勋视角 年9月15日...
    茶点故事阅读 34,518评论 2赞 309
  • 白月光启示录
    正文 我和宋清朗相恋三年，在试婚纱的时候发现自己被绿了。 大学时的朋友给我发了我未婚夫和他白月光在一起吃饭的照片。...
    茶点故事阅读 36,305评论 1赞 326
  • 活死人
    序言：一个原本活蹦乱跳的男人离奇死亡，死状恐怖，灵堂内的尸体忽然破棺而出，到底是诈尸还是另有隐情，我是刑警宁泽，带...
    沈念sama阅读 32,190评论 3赞 312
  • 日本核电站爆炸内幕
    正文 年R本政府宣布，位于F岛的核电站，受9级特大地震影响，放射性物质发生泄漏。R本人自食恶果不足惜，却给世界环境...
    茶点故事阅读 37,550评论 3赞 298
  • 男人毒药：我在死后第九天来索命
    文/蒙蒙 一、第九天 我趴在偏房一处隐蔽的房顶上张望。 院中可真热闹，春花似锦、人声如沸。这庄子的主人今日做“春日...
    开封第一讲书人阅读 28,880评论 0赞 17
  • 一桩弑父案，背后竟有这般阴谋
    文/苍兰香墨 我抬头看了看天上的太阳。三九已至，却和暖如春，着一层夹袄步出监牢的瞬间，已是汗流浃背。 一阵脚步声响...
    开封第一讲书人阅读 30,152评论 1赞 250
  • 情欲美人皮
    我被黑心中介骗来泰国打工， 没想到刚下飞机就差点儿被人妖公主榨干…… 1. 我叫王不留，地道东北人。 一个月前我还...
    沈念sama阅读 41,451评论 2赞 341
  • 代替公主和亲
    正文 我出身青楼，却偏偏与公主长得像，于是被迫代替她去往敌国和亲。 传闻我的和亲对象是个残疾皇子，可洞房花烛夜当晚...
    茶点故事阅读 40,637评论 2赞 335