测序数据比对到参考基因组

基因课FTP地址:ftp://http://gsx.genek.tv/2020-3-10%E7%9B%B4%E6%92%AD%E4%B8%80%E4%B8%AA%E5%AE%8C%E6%95%B4%E7%9A%84%E8%BD%AC%E5%BD%95%E7%BB%84%E9%A1%B9%E7%9B%AE/
听张旭东老师的课

服务器间数据拷贝

两台服务器间的数据拷贝 用 scp 用户名@服务器:文件路径

样本名称处理

  • sed -i 's/.fastq.gz//' xxx.file,类似rename的方式处理文件内容
  • awk 利用文件名生成样本信息表

构建参考基因组

hisat2-build genome.fasta genome
运行时间 —— 二十分钟以内

比对

hisat2 -x [ref-genome] -U [input_filename].fq.gz -S [output_file name].sam -p [threads] --new-summary --rna-strandness R

  • -U 写基因组名,不写后缀
  • PE 的输入文件-U项换为 -1 和 -2
  • 线程数建议一般设置为2-6即可
  • --new-summary 历史原因,使用tophat的旧日志文件格式则不加此项,用新格式日志文件则加此项
  • --rna-strandness 链特异性文库
    链特异性测序针对性解决的问题是:某些基因所在的正链与另一些基因的反链有交集,表达量定给谁?
    若不是链特异性测序,去掉此项
    若是链特异性测序,要问清楚是用的哪个技术, 大部分都用的是dUTP(90%) 如果是,
    单末端测序(SE) --rna-strandness参数设置为R
    PE 设置为RF
  • 目前绝大部分为链特异性测序
  • 非连特异性测序按照链特异性测序比对,有问题
    连特异性测序按照非链特异性测序比对,问题不大
  • 对于网上下载的测序数据有一些没有写明是链特异性测序还是非链特异性测序,解决方法:先假设是非连特异性测序,比对后在IGV上发现序列都是同一方向,则为链特异性文库
  • 比对DNA序列到基因组用bwa软件
  • 批量生成比对脚本,用awk实现
    vim的Ctrl+v也可以实现
    linux for循环也可以实现
  • 一般不需要设置错配率,默认就好。若比对后发现比对率特别低,则需要考虑。
  • 比对率一般至少70以上,比对率和 参考基因组测序组装质量、比对软件、测序品种与参考基因组物种亲缘关系 相关
  • 并行总线程可超过CPU数,超过即排队

比对结果查看

PE比对结果log文件

比对率97.44%

比对结果比对率统计与可视化

比对率结果在.log文件中
用软件MultiQC将多个log文件的结果统计

比对结果压缩排序

samtools sort -o xxx.bam xxx.sam

  • 这步比较耗内存,可以一个一个来
  • samtools view是只负责sam → bam的格式转换
  • bam文件构建好后,sam文件就可以删除了

对一个bam文件进行统计

samtools flagstat xxx.bam

  • 统计比对率(不同软件比对出来有差异)
  • 不同比对策略算出来的比对率略有不同,有primery的比对和secondary的比对,有区别
  • 转录组中建议以hisat2统计结果为准
    samtools的统计是通用的,没有对特定软件进行优化

构建bam index

samtools index xxx.bam

  • 构建index在转录组分析中除了IGV展示,没有其他用处
  • 在重测序中,bam文件构建index可用于变异检测等

IGV可视化

  • IGV官网 software.broadinstitute.org/software/igv/download
  • 用Java写的软件,跨平台(优点),易报错(劣)
  • 需要的文件
    ①导入基因组文件 genome.fasta
    ②基因注释文件genes.gtf
    ③sample.bam
    ④sample.bam.bai
  • 使用步骤
    • 建立基因组库
      Genomes → Creat .genome File...
    • 加载bam文件
      File → Load from file


      igv基因组建立

      gtf文件内容差不多长这样,只有transcript和exon
igv查看mapping结果
  • 桌面软件
  • 同一基因区同一碱基处,若大概一半一半的A/C概率,则为杂合,若极少量SNP可能是测序错误或RNA编辑
  • DNA测序鉴定突变:参考基因组为A,测序基因组绝大部分是A
  • 个体重测序更关心的是基因分型,不是只关心变异区域,表型=基因+环境,表型不完全有基因型决定。关心的是比例问题,概率问题

附加:转录本基因结构组装

  • hisat2对应软件为stringtie
  • 如果别人的基因组组装做的太差,可选择自己组装
  • 自己组装建议使用PASA流程,组装转录本

代码集中营

nohup hisat2-build xxx_genome.fasta xxx_genome 1>hisat2-build.log 2>&1 & # 标准输出与错误输出到同一文件
# 比对
# SE
hisat2 -x [ref-genome] -U [input_filename].fq.gz -S [output_file name].sam -p [threads] --new-summary --rna-strandness R 1>hisat2.log 2>&1
# PE
hisat2 -x [ref-genome] -1 [input_filename]_1.fq -2 [input_filename]_2.fq -S [output_file name].sam -p [threads] --new-summary --rna-strandness RF 1>hisat2.log 2>&1

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