2020转录组RNA-SEQ上游分析

文章目录

- 安装配置conda
- - 使用清华源下载sh脚本并安装
  - 设置镜像源
  - conda环境创建
  - 激活和退出环境
  - conda安装软件
- 质量评估 @ fastQC
- - fastq格式
  - FsatQC软件
  - multiqc综合所有qc
  - 结果解读
  - - 单一碱基占比
    - 单碱基测序质量在150bp上的分布情况
    - 测序质量的分布图
    - GC含量测定
    - 接头adapter统计
- 过滤reads @ trim_galore
- - 过滤标准
  - trim_galore
- 比对 @ hisat2
- - 参考基因组，gtf文件选择，不同基因组版本差异
  - 构建索引
  - - hisat的索引策略
    - hisat的比对策略
    - 索引的构建
  - hisat2比对
- sam格式转换 @ samtools
- - 头部区
  - 排序和转BAM
  - - 排序转bam
    - 索引bam文件
- 计数 @ featureCounts
- - 主要参数
  - 例子：
  - 结果处理

安装配置conda

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

# 使用清华源下载sh脚本
wget -c  /anaconda/miniconda/Miniconda3-latest-Linux-x86_64.sh

# 从官网下载最新版Miniconda3安装包，但速度较慢
wget -c /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 /anaconda/pkgs/free/
conda config --add channels /anaconda/pkgs/main/
conda config --add channels /anaconda/cloud/conda-forge/
conda config --add channels /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:
  - /anaconda/cloud/bioconda/
  - /anaconda/cloud/conda-forge/
  - /anaconda/pkgs/main/
  - /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格式描述：/s/8g-oUjiEhV4cGMJNuhmISQ
FastQ格式wiki：https://en.wikipedia.org/wiki/FASTQ_format
FastQ格式文献：/pmc/articles/PMC2847217/

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

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

1.第一行：必须以“@”开头，后面跟着唯一的序列ID标识符，然后跟着可选的序列描述内容，标识符与描述内容用空格分开；

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

FsatQC软件

FastQC质量评估软件官网：/projects/fastqc/
注意⚠️

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

常用参数：

# 常用参数
fastqc -o <> -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 <> *.

结果解读

单一碱基占比

unique reads 的占比，

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

也有这样的

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

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

测序质量的分布图

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

GC含量测定

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

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

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

接头adapter统计

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

过滤reads @ trim_galore

过滤标准

一般的过滤标准：

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

trim_galore

trim_galore使用手册：/FelixKrueger/TrimGalore/blob/master/Docs/Trim_Galore_User_Guide.md
阅读延伸：/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环境来运行
可以不用解压也可以运行

# 常用参数
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

批处理

ls * >
ls * >
paste   >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 &

批处理参考：/p/5a77b4156a5d
并行处理要求：/p/74a328f4d23a

比对 @ hisat2

hisat用户手册：/software/hisat2/
简书：有关hisat深入介绍，推荐：/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:///pub/ 。可以选择最新的 101版本（截止2020/08/25）。

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

不同基因组版本差异：

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

参考&资料：

/article/2017/07/2604/

bioinformatic论坛关于mask版本的讨论

一个介绍不是很清楚的不同基因组差异

构建索引

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

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

索引的构建

注意，如果要用到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..primary_assembly.fa homo_GRCh38. \
&

索引构建完毕如下：

hisat2比对

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

-x <path/to/base>             #索引的前缀,也可以添加路径，例如 ～/genome/homo_GRCh38.
-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 <> -x <homo_GRCh38_dna_primary> -1 <S_CTRL_1_val_1.>,<S_OE_1_val_1.> -2 <S_CTRL_2_val_2.>,<S_OE_2_val_2.> -S <./sam/>

批处理

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

# cat hista_align.sh
cat ~/lab/ZhangJinRui/analysis/trim_galore/ | while read id                                
do
    echo "hista2 :work_ID ${id} start at $(date)"
    echo "hisat2 -p 10 -q --known-splicesite-infile ~/lab/ZhangJinRui/analysis/geom/ -x homo_GRCh38_dna_primary -1 ~/lab/ZhangJinRui/analysis/trim_galore/${id}_1_val_1. -2 ~/lab/ZhangJinRui/analysis/trim_galore/${id}_2_val_2. -S ~/lab/ZhangJinRui/analysis/sam/${id}."
    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  > 
1

sam格式转换 @ samtools

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

头部区

@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格式信息：/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官方文档：/doc/
注意⚠️：

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

排序转bam

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

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

批处理

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

# cat sam_sort.sh 
cat |while read id  
do
    echo -e "*******WORK_ID ${id} START AT : \n$(date)"
    samtools sort -@ 10 -o ./sort_bam/${id}. ${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   #生成文件默认是原文件名 + .bai 也可以指定在原文件后面直接指定

批处理

ls *.bam > 

#cat bam_index.sh
cat |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

英文简明说明书：/featureCounts/
官方网站：/
英文用户手册下载：/subread-package/ （featurecounts说明在P31）
简书详细介绍（推荐）：/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信息解读：
/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] <>

     <>       #输入文件，sam/bam 支持多个文件输入
     -a <gtf>           #参考gtf注释文件，支持gz压缩
     -F <参考文件格式>    #默认GTF
     -T <int>           #线程数 1-32
     -J                 #对可变剪接计数
     -G <str>           #有-J设置时，用来提供参考基因组辅助寻找可变剪接
     -M                 #如果设置了-M则多重比对会被统计
     -o <>      #输出文件的名字，输出文件的内容是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  -p -t exon -g gene_id   *.

多个bam文件

featureCounts -t exon -g gene_id -a  -o

双端测序数据

featureCounts -p -t exon -g gene_id -a  -o  mapping_results_PE.bam

排除映射到不同染色体的fragment/read

featureCounts -p -C -t exon -g gene_id -a  -o  mapping_results_PE.bam

只计数那些双端都匹配的reads

featureCounts -p -B -t exon -g gene_id -a  -o  mapping_results_PE.bam

个人例子

#输出gene_id
nohup featureCounts -T 10 -a ~/lab/ZhangJinRui/analysis/geom/Homo_sapiens.GRCh38. -o  -p -t exon -g gene_id ~/lab/ZhangJinRui/analysis/sam/sort_bam/*. &

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

结果处理

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

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

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

因为我的笔记分布****、简书、知乎专栏等比较零散，管理起来比较麻烦，思考再三申请了一个微信公众号，会更加方便地发布更多有关生信息、统计方面内容，如果你觉得有需要可以尝试性和我缔结契约 ^ . ^ 。公众号如下：

秒客网

2020转录组RNA-SEQ上游分析

文章目录

安装配置conda

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

设置镜像源

conda环境创建

激活和退出环境

conda安装软件

质量评估 @ fastQC

fastq格式

FsatQC软件

multiqc综合所有qc

结果解读

单一碱基占比

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

测序质量的分布图

GC含量测定

接头adapter统计

过滤reads @ trim_galore

过滤标准

trim_galore

比对 @ hisat2

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

构建索引

hisat的索引策略

hisat的比对策略

索引的构建

hisat2比对

sam格式转换 @ samtools

头部区

排序和转BAM

排序转bam

索引bam文件

计数 @ featureCounts

主要参数

例子：

结果处理

相关文章