转录组练习(5)

这个部分主要是序列比对和reads读取,别人的帖子写的很全面,但有点复杂,需要兼顾来看。

//www.greatytc.com/p/ff585b72f04e

http://www.biotrainee.com/home.php?mod=space&uid=424&do=thread&view=me&from=space

总体来说,就是首先,在HIST网站下载人类的序列(我理解就是index),可以用终端下,也可以在HIST网站下,在网站下会快很多,然后放到终端解压。然后就是之前转换好的fastq格式(有1和2的那个),经过index比对,形成sam文件,详细过程看下面命令,每个编号转换为sam大概35分钟左右(根据大小和计算机水平),这个sam真心大,随便一个都十几个G。转换完sam之后,需要把它转换为bam格式,然后序列重整之后再转回sam.count格式。然后就可以reads了,经过reads之后就会变成count了,就可以导入R里面做差异基因等进一步分析了。


比对软件很多,首先大家去收集一下,因为我们是带大家入门,请统一用hisat2,并且搞懂它的用法。 直接去hisat2的主页下载index文件即可,然后把fastq格式的reads比对上去得到sam文件。 接着用samtools把它转为bam文件,并且排序(注意N和P两种排序区别)索引好,载入IGV,再截图几个基因看看! 顺便对bam文件进行简单QC,参考直播我的基因组系列。

下载 index

axel ftp://ftp.ccb.jhu.edu/pub/infphilo/hisat2/data/hg19.tar.gz

tar -zxvf hg19.tar.gz     #解压

ps.可以再HIST下,对应物种,老鼠就下老鼠的。。。

正式比对

命令行目录:Desktop/zhuanlu_files

hg19目录:Desktop/zhuanlu_files/hg19

RNA-Seq Data: Desktop/zhuanlu_files/RNA-Seq

注释:

-t 记录时间

-x hg19(index)文件路径

-1 -2 测序的两个fastq文件(单端测序用-U)

-S 比对结果输出路径

mkdir -p RNA-Seq/aligned

for i in`seq 56 58`

do

hisat2 -t -x hg19/genome -1 RNA-Seq/SRR35899${i}.sra_1.fastq.gz -2 RNA-Seq/SRR35899${i}.sra_2.fastq.gz -S RNA-Seq/aligned/SRR35899${i}.sam

done

上步就是把fastq格式转换为sam格式,并列序列比对

比对解读

总共比对了28856780条,7.4%一次都没有比对到,77.18%比对到一次,15.41%大于一次

没有匹配上的7.4%,不按照顺序再匹配,有4.3%匹配了

之后把剩余部分用单端数据进行比对,64.28%没有比对上,26.16比对了一次,9.56%大于一次

总共比对到95.45%

SAMtools

SAMtools 详解:

//www.greatytc.com/p/15f3499a6469

SAM 数据格式是目前高通量测序中存放比对数据的标准格式,但是SAM 文件都很大,非常占空间,所以需要转到bam文件,而用的就是SAMtools软件。

注释

-S 输入sam文件

-b 指定输出的文件为bam

排序默认是根据染色体的位置

for i in`seq 59 62`

do

samtools view -S SRR35899${i}.sam -b > SRR35899${i}.bam

samtools sort -n SRR35899${i}.bam -o SRR35899${i}_nsorted.bam

samtools view -h SRR35899${i}_nsorted.bam > SRR35899${i}_count.sam

done

上步就是先将sam 文件转成bam文件,bam 文件用 read name (-n)排序,再转回到sam文件,把这些转换好的文件放上层文件夹也行,新建文件夹也行,直接放上层文件夹会方便一些。

mkdir -p RNA-Seq/matrix/

htseq-count ~/SRR/SRR3589918_count.sam ~/gencode.v27lift37.annotation.gtf > 18.count

Genecode/gencode.v27lift37.annotation.gtf 这个是之前,在Genecode网站上下好的注释,老鼠的注释可以看看原文。把gtf这个文件放入Genecode这个文件夹,下面是原文的命令代码。

htseq-count RNA-Seq/aligned/SRR3589959_count.sam mouse/gencode.vM10.annotation.gtf

htseq-count RNA-Seq/aligned/SRR3589960_count.sam mouse/gencode.vM10.annotation.gtf

htseq-count RNA-Seq/aligned/SRR3589961_count.sam mouse/gencode.vM10.annotation.gtf

htseq-count RNA-Seq/aligned/SRR3589962_count.sam mouse/gencode.vM10.annotation.gtf

之后就可以得到count的文件,进行下一步处理。

最后编辑于
©著作权归作者所有,转载或内容合作请联系作者
  • 序言:七十年代末,一起剥皮案震惊了整个滨河市,随后出现的几起案子,更是在滨河造成了极大的恐慌,老刑警刘岩,带你破解...
    沈念sama阅读 211,884评论 6 492
  • 序言:滨河连续发生了三起死亡事件,死亡现场离奇诡异,居然都是意外死亡,警方通过查阅死者的电脑和手机,发现死者居然都...
    沈念sama阅读 90,347评论 3 385
  • 文/潘晓璐 我一进店门,熙熙楼的掌柜王于贵愁眉苦脸地迎上来,“玉大人,你说我怎么就摊上这事。” “怎么了?”我有些...
    开封第一讲书人阅读 157,435评论 0 348
  • 文/不坏的土叔 我叫张陵,是天一观的道长。 经常有香客问我,道长,这世上最难降的妖魔是什么? 我笑而不...
    开封第一讲书人阅读 56,509评论 1 284
  • 正文 为了忘掉前任,我火速办了婚礼,结果婚礼上,老公的妹妹穿的比我还像新娘。我一直安慰自己,他们只是感情好,可当我...
    茶点故事阅读 65,611评论 6 386
  • 文/花漫 我一把揭开白布。 她就那样静静地躺着,像睡着了一般。 火红的嫁衣衬着肌肤如雪。 梳的纹丝不乱的头发上,一...
    开封第一讲书人阅读 49,837评论 1 290
  • 那天,我揣着相机与录音,去河边找鬼。 笑死,一个胖子当着我的面吹牛,可吹牛的内容都是我干的。 我是一名探鬼主播,决...
    沈念sama阅读 38,987评论 3 408
  • 文/苍兰香墨 我猛地睁开眼,长吁一口气:“原来是场噩梦啊……” “哼!你这毒妇竟也来了?” 一声冷哼从身侧响起,我...
    开封第一讲书人阅读 37,730评论 0 267
  • 序言:老挝万荣一对情侣失踪,失踪者是张志新(化名)和其女友刘颖,没想到半个月后,有当地人在树林里发现了一具尸体,经...
    沈念sama阅读 44,194评论 1 303
  • 正文 独居荒郊野岭守林人离奇死亡,尸身上长有42处带血的脓包…… 初始之章·张勋 以下内容为张勋视角 年9月15日...
    茶点故事阅读 36,525评论 2 327
  • 正文 我和宋清朗相恋三年,在试婚纱的时候发现自己被绿了。 大学时的朋友给我发了我未婚夫和他白月光在一起吃饭的照片。...
    茶点故事阅读 38,664评论 1 340
  • 序言:一个原本活蹦乱跳的男人离奇死亡,死状恐怖,灵堂内的尸体忽然破棺而出,到底是诈尸还是另有隐情,我是刑警宁泽,带...
    沈念sama阅读 34,334评论 4 330
  • 正文 年R本政府宣布,位于F岛的核电站,受9级特大地震影响,放射性物质发生泄漏。R本人自食恶果不足惜,却给世界环境...
    茶点故事阅读 39,944评论 3 313
  • 文/蒙蒙 一、第九天 我趴在偏房一处隐蔽的房顶上张望。 院中可真热闹,春花似锦、人声如沸。这庄子的主人今日做“春日...
    开封第一讲书人阅读 30,764评论 0 21
  • 文/苍兰香墨 我抬头看了看天上的太阳。三九已至,却和暖如春,着一层夹袄步出监牢的瞬间,已是汗流浃背。 一阵脚步声响...
    开封第一讲书人阅读 31,997评论 1 266
  • 我被黑心中介骗来泰国打工, 没想到刚下飞机就差点儿被人妖公主榨干…… 1. 我叫王不留,地道东北人。 一个月前我还...
    沈念sama阅读 46,389评论 2 360
  • 正文 我出身青楼,却偏偏与公主长得像,于是被迫代替她去往敌国和亲。 传闻我的和亲对象是个残疾皇子,可洞房花烛夜当晚...
    茶点故事阅读 43,554评论 2 349

推荐阅读更多精彩内容