转录组那些事儿 Part II

刘小泽写于18.8.21
争取用三部分整理完
Part I是前情提要，学习一些背景知识；
Part II是实战，从文章拿到数据开始，只要是服务器能完成的工作都完成；
Part III 将是R处理部分，也是最核心的部分，我会先学习Bioconductor，抓紧时间吧

. 阅读文献

文章梗概

运动神经元存活蛋白(Survival Motor Neuron, SMN)是人体内的泛表达蛋白，它的会导致脊髓性肌萎速症(Spinal muscular atrophy, SMA)，SMA是由SMN1基因突变引起的，但是这个基因是广泛表达的基因，为什么运动神经元(Motor Neurons, MNs)偏偏是最易受影响的细胞类型之一？实验通过对照组和SMA患者诱导多能干细胞（IPSC）从而产生纯化的运动神经元，通过固定、抗体标记，然后进行了RNA测序研究。在患者运动神经元中发现了SMA特异性的变化，其中包括内质网(ER)应激通路的过度激活。功能研究表明，抑制内质网的应激反应能提高运动神经元的存活率。在患病小鼠中，使用ER应激抑制剂穿过血脑屏障可以保护脊髓运动神经元。因此，选择性激活内质网应激反应，导致了SMA患病个体的运动神经元死亡。

实验流程
GEO数据

开始实战

万事之源，软件为先

#配置conda
conda create -n rna-seq python=2 samtools fastp sra-tools hisat2 rseqc -y
conda install -c bcbio htseq -y
conda install numpy pysam -y
CONBIN=/home/biosoft/miniconda3/envs/rna-seq/bin

配置好工作路径

RNA=/home/project/rna-seq
mkdir -p $RNA/{raw_data,clean_data,ref/{genome,gtf,index},align,stats,matrix}
RAW=$RNA/raw_data
CLEAN=$RNA/clean-data
GENOME=$RNA/ref/genome
GTF=$RNA/ref/gtf
INDEX=$RNA/ref/index
ALIGN=$RNA/aign
MATRIX=$RNA/matrix
STATS=$RNA/stats

下载测试数据

GEO数据库中的编号是：GSE69175

关于数据：

NGS测序数据一般会上传到SRA数据库里面，但是怎么从GEO数据库中找到SRA原始数据的下载地址？【GEO数据库地址：https://www.ncbi.nlm.nih.gov/geo/】

GEO数据库下载

具体的层级关系是：SRP(项目)—>SRS(样本)—>SRX(数据产生)—>SRR(数据本身)

SRR数据库地址：https://ftp-private.ncbi.nlm.nih.gov/sra/sra-instant/reads/ByRun/sra/SRR/

##################################################################
################# Part 1: Data download #########################
##################################################################
# SRR2038215-SRR2038216: 对照
# SRR2038217-SRR2038218: SMN
cd $RAW
for i in `seq 15 18`;do 
ascp -v -i ~/.aspera/connect/etc/asperaweb_id_dsa.openssh \
-k 1 -T -l200m anonftp@ftp-private.ncbi.nlm.nih.gov:/sra/sra-instant/reads/ByRun/sra/SRR/SRR358/SRR35899${i}/SRR35899${i}.sra ./ #使用ascp下载，嗖嗖嗖～
done
fastq-dump --gzip --split-3 *.sra #sra转为fastq
#以下3行原版来自生信技能树Jimmy
find $RAW -name "*.gz" | grep 1.fastq.gz >1
find $RAW -name "*.gz" | grep 2.fastq.gz >2
paste 1 2 > raw_conf #将read1和read2各自的合集再整合到一起，形成一个数据配置文件
cp raw_conf $CLEAN

质控过滤

这里使用的是fastp，同时融合了质控、过滤的模式；
同时也可以使用fastqc + trimmomatic/ trim_galore进行

##################################################################
################## Part 2: QC & trim ############################
##################################################################
source activate rna-seq
cd $CLEAN
#下面四行原版来自生信技能树Jimmy
cat raw_conf | while read id;do
line=(${id})
fq1=${line[0]}
fq2=${line[1]}

fastp -i $fq1 -I $fq2 -o out.$(basename $fq1) -O out.$(basename $fq2)-w 16
done
source deactivate

下载参考基因组和注释文件

人类参考基因组选择UCSC数据库；注释文件选择GENCODE,https://www.gencodegenes.org/

下载好基因组，需要构建基因组索引

如果自己的项目做非模式物种，可以用二代三代测序数据组装成参考转录组，例如trinity组装。如果做的物种有参考基因组，就方便一些，可以直接从相关数据库中下载参考基因组

##################################################################
############### Part 3: prepare ref & index ######################
##################################################################

##Download hg19 
cd $GENOME
for i in $(seq 1 22) X Y M;do
wget http://hgdownload.cse.ucsc.edu/goldenPath/hg19/chromosomes/chr${i}.fa.gz;
done
gunzip *.gz
for i in $(seq 1 22) X Y M;do
cat chr${i}.fa >> hg19.fa
done
rm chr*
##或者也可以直接下载成品
#wget http://10.10.0.8/hgdownload.soe.ucsc.edu/goldenPath/hg19/bigZips/hg19.fa.gz

##Download GRCh38 https://www.gencodegenes.org/releases/current.html
cd $GTF
wget ftp://ftp.ebi.ac.uk/pub/databases/gencode/Gencode_human/release_28/GRCh37_mapping/gencode.v28lift37.annotation.gtf.gz
gunzip *.gz
## hisat2 index 
cd $INDEX
hisat2-build -p 20 $GENOME/hg19.fa hg19 #-p为线程数

比对

##################################################################
######## Part 4: align & sort & index human samples:56-58 ########
##################################################################
source activate rna-seq
cd $ALIGN
for i in `seq 15 18`;do
hisat2 -t -p 20 -x $INDEX/hg19 \
-1 $CLEAN/out.SRR20382${i}_1.fastq.gz \
-2 $CLEAN/out.SRR20382${i}_2.fastq.gz \
-S SRR20382${i}.sam
samtools view -Sb SRR20382${i}.sam > SRR20382${i}.bam
samtools sort -@ 20 -o SRR20382${i}.sorted.bam SRR20382${i}.bam
samtools index SRR20382${i}.sorted.bam
rm *.sam
done
source deactivate
#关于排序：默认按染色体位置排序；-n根据read名排序；-t根据tag排序

使用了samtools的三件套：转换(view)、排序(sort)、建索引(index)
转换: -S指输入文件格式（不加-S默认输入是bam），-b指定输出文件（默认输出sam）【如果要bam转sam，-h设置输出sam时带上头注释信息】
排序： 对bam排序，-@ 线程数， -o输出文件
索引： 结果会产生.bai文件【必须排序后才能建索引～就像体育课必须从高到矮排好以后再报数】

基本信息统计

##################################################################
################ Part 5: basic statistics #######################
##################################################################
source activate rna-seq
cd $STATS
for i in `seq 15 18`;do
samtools flagstat $ALIGN/SRR20382${i}.sorted.bam > SRR20382${i}.sorted.flag
done
#如果想根据flag提取特定区域，比如想查看1号染色体100-10000的区域的信息
#samtools view -b -f 0x10 $ALIGN/SRR20382${i}.sorted.bam chr1:100-10000 > ${i}.flag.bam
#samtools flagstat ${i}.flag.bam

#使用RSeqQC统计
#先用bam_stat.py对bam文件统计，看下比对率
#bam_stat.py -i $ALIGN/SRR20382${i}.sorted.bam > SRR20382${i}.bam.stat
source deactivate
基本信息统计

reads计数

基于基因组水平，可以用Htseq-count和featureCounts

  1. Htseq-count：它是用python写的一个脚本，conda install -c bcbio htseq -y安装好以后可以直接拿来用

    需要比对文件sam/bam格式；需要基因组注释文件gff/gtf格式。这两个文件的染色体名称必须相同，因为需要将比对位置和特征信息相匹配就要借助染色体名。
    它的工作原理是：先通过bam文件找到reads比对上的外显子，然后去gtf文件中将外显子的基因ID进行统计，得到的结果就是未经标准化的表达量

    source activate rna-seq
cd $MATRIX
for i in `seq 15 18`;do
$CONBIN/htseq-count -s no -r name -f bam $ALIGN/SRR20382${i}.sorted.bam \
$GTF/gencode.v28lift37.annotation.gtf \
>SRR20382${i}.count 2>SRR20382${i}.log
done
source deactivate

    htseq-count的-s参数是strand意思，默认使用链特异性建库方式，也就是计算过程中要求read1只能和双端测序的其中一条比较，read2只能和另一条比较；我们一般设置no即可，表示每条read都和正负链比较一次；
    -r参数 两个选择name和pos ，数据根据位置或者read名称排序，默认name；
    -f参数 输入文件格式bam/sam；

    -m参数一般也就是默认union模式

    结果每个样本输出一个count文件，其中包含了基因名和reads计数；另外，如果你看到文件倒数几行，会发现还有几行带文字的
    htseq-count计数

    no_feature：比对区域与任何基因都没有重叠
    ambiguous：比对区域与多个基因都发生重叠
    too_low_aQual：比对质量低于设定阈值（默认是10）
    not_aligned：无法比对上
    alignment_not_unique：比对位置不唯一

  2. featureCounts：隶属于subread套件【相比于htseq更快】

    source activate rna-seq
cd $MATRIX
begin=$(date +%s)
for i in `seq 15 18`;do
featureCounts -T 20 -p -t exon -g gene_id -a $GTF/gencode.v28lift37.annotation.gtf  -o SRR20382${i}.feature.count $ALIGN/SRR20382${i}.bam
done
tim=$(echo "$(date +%s)-$begin" | bc)
printf "time of featureCounts for 4 samples: %.4f seconds" $tim
source deactivate
#运行结果显示：
#time of featureCounts for 4 samples: 366.5310 seconds

    -T参数设置线程，
    -p参数表示针对双端测序数据，
    -t 参数默认将外显子认定为基因组的一个feature(搜索特征)，
    -g参数 默认按照基因名来计数，
    -a参数 输入注释文件，
    -o参数 输出文件

    结果有两个，一个count文件，一个summary文件


    featureCounts文件

对两个软件的结果进行合并

##合并htdeq生成的count文件到matrix.count
cd $MATRIX/htseq
perl -lne 'if ($ARGV=~/(.*).count/){print "$1\t$_"}' *.count | grep -v "_" >matrix.count
##合并featureCounts生成的count文件到matrix_2.count
cd $MATRIX/featureCounts
for i in `seq 15 18`;do
cut -f 1,7 SRR20382${i}.feature.count | grep -v "^#" > SRR20382${i}.matrix
sed -i '1d' SRR20382${i}.matrix
done
perl -lne 'if ($ARGV=~/(.*).matrix/){print "$1\t$_"}' *.matrix >matrix_2.count

统计一下两个软件的计数之和

#统计featureCounts
perl -alne '$sum += $F[2]; END {print $sum}' matrix.count
#结果是5882943
#统计htseq-count,结果是786338

欢迎关注我们的公众号～_～　　
我们是两个农转生信的小硕，打造生信星球，想让它成为一个不拽术语、通俗易懂的生信知识平台。需要帮助或提出意见请后台留言或发送邮件到Bioplanet520@outlook.com

Welcome to our bioinfoplanet!

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

推荐阅读更多精彩内容