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

前期文章

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

3 Alignment

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

首先学一下hisat2的用法
hisat2 --help  ## 主要要学会使用--help等参数调用帮助文档(软件说明书)
Usage:
hisat2 [options]* -x <ht2-idx> {-1 <m1> -2 <m2> | 
-U <r> | --sra-acc <SRA accession number>} [-S <sam>]

-p/--threads <int> number of alignment threads to launch (1)
用法一,一步一个脚印
## 用“指鹿为马“的方法,将会重复出现的路径等缩略
wkd=/Users/Desktop/project

## 比对
cd $wkd/clean 

ls *gz|cut -d"_" -f 1,2 |sort -u |\
while read id; 
do 
   ls -lh ${id}R1_val_1.fq.gz ${id}R2_val_2.fq.gz 
   hisat2 -p 4 -x \
          /Users/Desktop/project/reference/hg38/hisat2 \ 
          -1 ${id}1_val_1.fq.gz \
          -2 ${id}2_val_2.fq.gz -S ${id}.hisat.sam 
done 

## 转为二进制文件sam to bam
ls *.bam| while read id; 
do 
(samtools sort -@ 4 -o $(basename ${id}".sam").bam ${id});
done
rm *.sam  # 移除sam文件
用法二:一步到位
vim align.sh ## 创建一个脚本,写入以下内容

#!/bin/bash
set -e
set -u
set -o pipefail

# source activate rna

# set PATH
HISAT2=/Users/Desktop/project/reference/hg38/hisat2

pwd

ls *gz | cut -d"_" -f 1,2 | sort -u |\
while read id
do
     echo "processing  ${id}_R1_val_1.fq.gz ${id}_R2_val_2.fq.gz"
     hisat2 -p 4 -x \
            $HISAT2/hg38.ht2 \ 
            -1 ${id}_R1_val_1.fq.gz \
            -2 ${id}_R2_val_2.fq.gz |\
            samtools sort -@ 4 > ${id}.hisat.bam


done
# source deactivate

运行脚本

nohup bash align.sh &

查看结果

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

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

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

小结:本章所需知识背景

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

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

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

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

问题集锦

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

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


    QC报错?

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

for id in {xx1,xx2,xxx3};
do 
echo "Processin sample ${id}"
hisat2 -p 5 -x /data1/reference/index/hisat2/hg19/hg19 \
-1 /data1/projects/clean/${id}_R1_trimmed.fq.gz \
-2 /data1projects/clean/${id}_R2_trimmed.fq.gz \
-S /data1/projects/align1/${id}.hisat.sam
done 
  1. 为什么比对率这么低?


    image.png

    解答:


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


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

感谢洲更

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