专栏名称: 生信媛
生信媛,从1人分享,到8人同行。坚持分享生信入门方法与课程,持续记录生信相关的分析pipeline, python和R在生物信息学中的利用。内容涵盖服务器使用、基因组转录组分析以及群体遗传。
目录
相关文章推荐
51好读  ›  专栏  ›  生信媛

一个超简单的转录组项目全过程--iMac+RNA-Seq(三)Alignment 比对

生信媛  · 公众号  · 生物  · 2019-07-05 17:06

正文

请到「今天看啥」查看全文



前期文章

一个超简单的转录组项目全过程--iMac+RNA-Seq(一)

一个超简单的转录组项目全过程--iMac+RNA-Seq(二)QC

3 Alignment

有很多软件都可以比对到参考基因组,hisat2,bowtie2,STAR等等(手上这台iMac不够内存跑STAR)

首先学一下hisat2的用法
  1. hisat2 --help ## 主要要学会使用--help等参数调用帮助文档(软件说明书)

  2. Usage:

  3. hisat2 [options]* -x <ht2-idx> {-1 -2 |

  4. -U | --sra-acc <SRA accession number>} [-S ]


  5. -p/--threads number of alignment threads to launch (1)

用法一,一步一个脚印
  1. ## 用“指鹿为马“的方法,将会重复出现的路径等缩略

  2. wkd=/Users/Desktop/project


  3. ## 比对

  4. cd $wkd/clean


  5. ls *gz|cut -d"_" -f 1,2 |sort -u |\

  6. while read id;

  7. do

  8. ls -lh ${id}R1_val_1.fq.gz ${id}R2_val_2.fq.gz

  9. hisat2 -p 4 -x \

  10. /Users/Desktop/project/reference/hg38/hisat2 \

  11. -1 ${id}1_val_1.fq.gz \

  12. -2 $ {id}2_val_2.fq.gz -S ${id}.hisat.sam

  13. done


  14. ## 转为二进制文件sam to bam

  15. ls *.bam| while read id;

  16. do

  17. (samtools sort -@ 4 -o $(basename ${id}".sam").bam ${id});

  18. done

  19. rm *.sam # 移除sam文件

用法二:一步到位
  1. vim align.sh ## 创建一个脚本,写入以下内容


  2. #!/bin/bash

  3. set -e

  4. set -u

  5. set - o pipefail


  6. # source activate rna


  7. # set PATH

  8. HISAT2=/Users/Desktop/project/reference/hg38/hisat2


  9. pwd


  10. ls *gz | cut -d"_" -f 1,2 | sort -u |\

  11. while read id

  12. do

  13. echo "processing ${id}_R1_val_1.fq.gz ${id}_R2_val_2.fq.gz"

  14. hisat2 -p 4 -x \

  15. $HISAT2/hg38.ht2 \

  16. -1 ${id}_R1_val_1.fq.gz \

  17. -2 ${id}_R2_val_2. fq.gz |\

  18. samtools sort -@ 4 > ${id}.hisat.bam



  19. done

  20. # source deactivate

运行脚本

nohup bash align.sh &

查看结果

请关注比对报告中的aligned concordantly exactly 1 time这一项内容,至少要80 %以上。

小彩蛋 :洲更老师的锦囊 jobs and disown的使用

当我忘记nohup脚本的时候,但我又要离开实验室了,这时候使用disown帮助把跑到一半的任务丢到后台,惊呆了!神仙走位!

小结:本章所需知识背景

(1)什么是比对,比对前后数据的变化;

(2)fasta,fastq文件的规则以及包含的信息;

(3)比对结果的查看,比对率的由来。

以上知识点均可以在网上查到,建议在进行着一部分项目的时候,一定要补充上最基本的知识点。

问题集锦

收到了几个小问题,不过我几乎都无法解答,但我可以找到答案~

我的怎么QC怎么回事?跑着跑着跑着就停了呀。

QC报错?

问了代码才知道问题可能出在input上,过滤结束后的样本名是val.fq.gz ,出现这样结果的原因,应该是过滤出错了。

  1. for id in {xx1,xx2,xxx3};

  2. do

  3. echo "Processin sample ${id}"

  4. hisat2 -p 5 -x /data1/reference/index/hisat2/hg19/hg19 \

  5. -1 /data1/projects/clean/${id}_R1_trimmed.fq.gz \

  6. -2 /data1projects/clean/${id}_R2_trimmed.fq.gz \

  7. -S /data1/projects/align1/${id}.hisat.sam

  8. done


为什么比对率这么低?

解答:

比对率的计算好像不对? 怎么加起来都不等于96.7%呀?这样的问题,一般我都解决不了,我会直接放弃,然后寻求帮助。

解答: 具体请看简书:2018-12-04 Hisat2 map 结果与 samtools flagstat 结果不一致








请到「今天看啥」查看全文