转录组分析流程：比对（有参）及统计Counts矩阵

1. 质控

fastqc *
multiqc *
trimmomatic_run.sh #去掉前9个碱基

trimmomatic_run.sh

#! /bin/bash
#history:
#   Gossie V1 20190101#用trimmomatic对进行测序数据质量控制
#输入文件为包含样本名称的单列文件#线程
t=1
# 模式PE（双端）或者SE（单端）
mode='PE'#以下参数适情况修改、添加。
#ILLUMINACLIP模式
#接头文件，fasta格式。注意：不同的命名格式有差别。
fastaWithAdapterEtc='/Trimmomatic-0.38/adapters/TruSeq3-PE-2.fa'
#第一步搜索时允许搞错错配碱基个数，通常为1-2.
seed_mismatches=2
#双端测序的palindrome模式下，连上接头的两条reads的比对阈值
palindrome_clip_threshold=30
#切除的接头序列的最低比对分值，通常为7-15.
simple_clip_threshold=10
#仅对palindrome模式有效，双端测序palindrome模式下，可以切除的接头的最短长度。默认是8。但，可以设置为1.
minAdapterlength=8
#仅对palindrome模式有效。默认为false。
keepBothReads=true# SLIDINGWINDOW:滑窗模式。从5‘端开始滑动，切除平均质量小于阈值的片段。
#窗口大小
windowSize=0
#质量
requiredQuality=0#LEADING:切除起位置小于阈值的碱基
lead_quality=0
#TRAILING:切除末尾位置小于阈值的碱基
trail_quality=0
#CROP:从末尾切除read部分碱基，使reads达到指定长度.
crop_length=150
#HEADCROP:切除起始位置特定长度的碱基
headcrop=9
#MINLEN:舍弃小于给定长度的
minlen=50
#AVGQUAL：舍弃平均质量低于给定水平的read
avgqual=20while read id
do
nohup trimmomatic PE -threads $t ${id}_1.fastq.gz ${id}_2.fastq.gz -baseout ${id}.fastq.gz HEADCROP:$headcrop ILLUMINACLIP:$fastaWithAdapterEtc:$seed_mismatches:$palindrome_clip_threshold:$simple_clip_threshold:$minAdapterlength:$keepBothReads MINLEN:$minlen &
done < $1

2. 比对

比对软件：hisat2

gffread annotation.gff3 -T -o annotation.gtf #hisat2需要GTF格式的注释文件，注意看一下有没有外显子信息(exon)，有的gff3可能不提供，会影响比对reads的统计
hisat2_extract_splice_sites.py .UTX.gtf >UTX.gene.ss #剪接位点文件hisat2_extract_exons.py UTX.gene.gtf >UTX.gene.exon #外显子位点文件sh hisat2.sh #samtools version 1.21

hisat2.sh
需要 1.9以上版本的samtools用于排序和转换bam

#! /bin/bash
# 输入文件：sra id
#index文件
index=/hisat2_index/UTX
while read id
do
#--dta 输出比对情况
nohup hisat2 -p 2 --dta -x $index -1 ${id}_1P.fastq.gz -2 ${id}_2P.fastq.gz -S $id.hisat2.sam >$id.hisat2.log 2>&1 && samtools sort -@ 10 -o $id.sorted.bam $id.hisat2.sam && rm $id.hisat2.sam &
done < $1

3. 计算count矩阵

软件：featureCounts


nohup featureCounts -p -T 64 -t exon -g gene_id -a UTX.gene.gtf -o counts.featurecounts.txt *.sorted.bam >featurecounts.log 2>&1 &
# -p：双端；单端不用加参数
# -T：线程，1~64
# -t：feature type
# -g：meta-feature 类型，若一个reads或fragment比对到多个feature（如，exon),如果这些exon属于一个meta-feature（如，gene)，则只计数一次
# -a：注释文件，最好用GTF；因为GFF3的ID注释包含gene、mRNA、exon、UTR等
# -o：输出文件
# 可对多个bam同时统计

转录组分析流程：比对（有参）及统计Counts矩阵相关推荐

转录组分析流程｜数据处理与De novo组装（一）
title: 转录组分析流程|数据处理与De novo组装(一) tags: - 转录组组装 - 教程 - 软件 - Trinity - Rcorrector - Trimmomatic catego ...
空间转录组分析流程（使用Seurat对空间数据集进行分析）
空间转录组分析流程(使用Seurat对空间数据集进行分析) 因为每次打开这个网页都非常慢,所以我讲这个网页进行一个翻译,方便学习. 使用Seurat对空间数据集进行分析,可视化和集成 1.介绍本教程 ...
生物信息学之rnaseq转录组分析流程--转换文件中的ensemble id到gene名
生物信息学之rnaseq转录组分析--转换文件中的ensemble id到gene名如何解决转录组分析中count之后遇到ensemble id的问题一个将ensemble id转换成gene名的 ...
转录组分析流程：表达差异分析之edgeR
edgeR edgeR是非常经典的转录组表达差异分析软件. 样本量:72个转录组样本 library(edgeR) library(HTSFilter)fc <- read.table('cou ...
转录组分析流程｜TransDecoder预测转录本的开放阅读框（二）
使用TransDecoder预测CDS TransDecoder按照其官网的说明,主要用于识别转录本序列中的潜在的编码区域,也就是预测CDS.转录本可以由RNA-Seq数据通过Trinity组装来的, ...
转录组分析流程｜基于salmon转录组批量定量流程（三）
TransDecoder那一步最终得到了*.cds序列,之后就需要用到salmon进行下面操作 salmon介绍 Salmon是不基于比对计数而直接对基因进行定量的工具,适用于转录组.宏基因组等的分析 ...
高级转录组分析和R语言数据可视化第十三期（线上线下同时开课）
" 福利公告:为了响应学员的学习需求,经过易生信培训团队的讨论筹备,现决定安排扩增子16S分析.宏基因组.Python课程线上直播课.报名参加线上直播课的老师可在1年内选择参加同课程的一次线 ...
本周开课 | 第 17 期高级转录组分析和R数据可视化火热报名中！！！
福利公告:为了响应学员的学习需求,经过易生信培训团队的讨论筹备,现安排<高级转录组分析和R数据可视化>于2022年4月15-17 线上/线下课程 (线上课是通过腾讯会议实时直播线下课,实时 ...
第16期高级转录组分析和R数据可视化培训(2022年1月)
福利公告:为了响应学员的学习需求,经过易生信培训团队的讨论筹备,现安排<高级转录组分析和R数据可视化>于2022年1月14-16 线上/线下课程 (线上课是通过腾讯会议实时直播线下课,实时 ...

转录组分析流程：比对（有参）及统计Counts矩阵

1. 质控

2. 比对

3. 计算count矩阵

转录组分析流程：比对（有参）及统计Counts矩阵相关推荐

最新文章

热门文章