跟着大神学生信 Jimmy RNA-Seq

2019-10-03 本文已影响0人致知_5974

暴躁版

制约我们的是热情吗？是知识吗？

不，是网速啊，阿西吧！

1 获取SRA号——自己读文献，去pubmend搜

2 下载数据

cat SRR_Acc_List.txt | while read id; do (prefetch  ${id} );done
别试了，龟速，下不下来，请愉快的打上一个&然后洗洗睡吧
jimmy贴心的给了一个kill命令的代码
ps -ef | grep prefetch | awk '{print $2}' | while read id; do kill ${id}; done
呵呵哒

3 质控

sra转fastq

for i in $wkd/*sra
do
        echo $i
        fastq-dump --split-3 --skip-technical --clip --gzip $i  ## 批量转换
done

核心：fastq-dump
--split-3 确保一个双端测序的样本被拆分成两个fastq文件。
--skip-technical Dump only biological reads
-W|--clip Remove adapter sequences from reads
--gzip Compress output using gzip: deprecated, not
recommended （哈哈，为啥不推荐呢）
fastq转fastqc

ls *gz | xargs fastqc -t 4
multiqc ./

-t 线程数（double哈哈）
还是好慢，rna-seq的正确打开方法是，找本书，一边seq一边看
multiqc一下，就得到了multiqc的报告，下下来看一下，打不开，loading report ，FU...K!马景涛咆哮版！

4.过滤

过滤质量差的reads和去除接头

ls ~/sra/*_1.fastq.gz >1
 ls ~/sra/*_2.fastq.gz >2
paste 1 2 >config

cat 一下config,长这样

/trainee2/hz10/sra/SRR1039510_1.fastq.gz    /trainee2/hz10/sra/SRR1039510_2.fastq.gz
/trainee2/hz10/sra/SRR1039511_1.fastq.gz    /trainee2/hz10/sra/SRR1039511_2.fastq.gz
/trainee2/hz10/sra/SRR1039512_1.fastq.gz    /trainee2/hz10/sra/SRR1039512_2.fastq.gz

trim_galore ,用于去除低质量和接头数据
改一下大神的脚本，换成自己的路径，长这样

conda activate rna2
bin_trim_galore=trim_galore
dir='/trainee2/hz10/sra/clean'
 cat $1 |while read id
do
       arr=(${id})
       fq1=${arr[0]}
       fq2=${arr[1]}
  $bin_trim_galore -q 25 --phred33 --length 36 --stringency 3 --paired -o $dir $fq1 $fq2
done
conda deactivate rna2

然后bash一下，好慢，我忘了打&，然后掉线了，我恨！我需要把之前trim好的删掉再重新bash吗？还是它会很智能的自己跳过去呢？它好像跳过去了，nice~
太慢了，回家吃饭了。。。。。
好了，如下。

├── [ 123]  1
├── [ 123]  2
├── [ 246]  config
├── [ 284]  qc.sh
├── [2.9K]  SRR1039510_1.fastq.gz_trimming_report.txt
├── [1.2G]  SRR1039510_1_val_1.fq.gz
├── [3.1K]  SRR1039510_2.fastq.gz_trimming_report.txt
├── [1.2G]  SRR1039510_2_val_2.fq.gz
├── [2.9K]  SRR1039511_1.fastq.gz_trimming_report.txt
├── [1.2G]  SRR1039511_1_val_1.fq.gz
├── [3.1K]  SRR1039511_2.fastq.gz_trimming_report.txt
├── [1.2G]  SRR1039511_2_val_2.fq.gz
├── [2.9K]  SRR1039512_1.fastq.gz_trimming_report.txt
├── [1.5G]  SRR1039512_1_val_1.fq.gz
├── [3.1K]  SRR1039512_2.fastq.gz_trimming_report.txt
└── [1.5G]  SRR1039512_2_val_2.fq.gz

--phred33 Instructs Cutadapt to use ASCII+33 quality scores as Phred scores (Sanger/Illumina 1.9+ encoding) for quality trimming.
--stringency Overlap with adapter sequence required to trim a sequence. 剪切掉的overlap的序列
***--length <INT> **** Discard reads that became shorter than length INT because of either quality or adapter trimming. A value of '0' effectively disables this behaviour. Default: 20 bp.
--paired 双端测序并且需要双端测序的两条序列都比length的设定值要长，可以再不影响fastq格式的情况下，过滤掉过短的序列。

4. 比对

使用hisat2
构建hisat2基因组索引,来自生信技能树论坛，传说很慢

# 构建目录
mkdir /mnt/d/reference/index
mkdir /mnt/d/reference/index/hisat
#下载
wget -P /mnt/d/reference/index/hisat [url=ftp://ftp.ccb.jhu.edu/pub/infphilo/hisat2/data/hg19.tar.gz]ftp://ftp.ccb.jhu.edu/pub/infphilo/hisat2/data/hg19.tar.gz[/url]
#解压，-C指定解压目录
tar -zxvf /mnt/d/reference/index/hisat/hg19.tar.gz -C /mnt/d/reference/index/hisat
# 移除下载的压缩包
rm /mnt/d/reference/index/hisat/*.tar.gz
#查看解压文件
ll /mnt/d/reference/index/hisat/hg19

Tips_1:不用单独给下载的hg19的index文件再分别构建目录，解压的时候会自动创建hg19目录。
Tips_2: 解压的文件中，包含genome.*.ht2的8个文件和一个shell脚本。
hisat2用法

hisat2 [options]* -x <hisat2-idx>{-1 <m1> -2 <m2> | -U <r> } [-S <hit>]

-x 指定基因组索引
-1 指定第一个fastq文件:双端测序结果的第一个文件。若有多组数据，使用逗号将文件分隔。Reads的长度可以不一致。
-2 指定第二个fastq文件:双端测序结果的第二个文件。若有多组数据，使用逗号将文件分隔，并且文件顺序要和-1参数对应。Reads的长度可以不一致。
-S 指定输出的SAM文件。
改了老师的批量运行代码

cd $wkd/clean 
ls *gz|cut -d"_" -f 1 |sort -u |while read id;do
ls -lh ${id}_1_val_1.fq.gz   ${id}_2_val_2.fq.gz 
hisat2 -p 10 -x ~/sra/index/hisat/hg38/genome -1 ${id}_1_val_1.fq.gz   -2 ${id}_2_val_2.fq.gz  -S ${id}.hisat.sam
done

运行了一下，运行完第一个不错，比对98%，第二个卡住了，报错，google了一下原因，sam文件太大，内存不够，服务器爆掉了。。。。我太难了。。。。
然后就一个一个做，运行完了转sam，或者直接转sam。

hisat2 -p 1 -x ~/sra/index/hisat2/hg38/genome  -1 SRR1039511_1_val_1.fq.gz  -2 SRR1039511_2_val_2.fq.gz  | samtools sort -@ 1 -o ~/sra/clean/SRR1039511.sort.bam - 
bam 转 sam
samtools sort -O bam -@ 5  -o SRR1039510.hisat.bam SRR1039510.hisat.sam

-O 输出文件格式
-@ 线程数
最后是要操作的文件名，上面用管道符的是最后的-
我想买个服务器，我太难了。。。。
终于结束了，为bam文件建立索引

ls *.bam |xargs -i samtools index {}
ls *.bam |while read id ;do ( samtools flagstat -@ 1 $id >  $(basename ${id} ".bam").flagstat  );done

samtools view SAM转换为二进制对应的BAM格式
samtools sort 比对排序
samtools index 构建索引。对排序文件进行索引之后，有利于快速提取基因组重叠区域的比对结果
samtools flagstats 给出BAM文件的比对结果

samtools常用命令详解

SRR1039512.sort.bam.bam.flagstat文件可以直接cat查看比对结果

5.获得表达矩阵

gtf="/teach/database/gtf/gencode.v29.annotation.gtf.gz" 
featureCounts -T 5 -p -t exon -g gene_id  -a $gtf -o  all.id.txt  ~/sra/clean/*.bam  1>counts.id.log 2>&1 &

featurecounts 据说最大的优点是快
-a 输入GTF/GFF基因组注释文件
-p 这个参数是针对paired-end数据。Check validity of paired-end distance when counting read pairs. Use -d and -D to set thresholds.
-F 指定-a注释文件的格式，默认是GTF
-g 从注释文件中提取Meta-features信息用于read count，默认是gene_id
-t 跟-g一样的意思，其是默认将exon作为一个feature
-o 输出文件
-T 多线程数
然后我们就得到表达矩阵 all.id.txt 啦，就可以用R进行下游分析啦啦啦啦（啦啦啦你妹），limma,DESeq2，想用什么用什么，哪里不会点哪里~
ps，我觉得jimmy大佬有一个地方写错了，align下面没有bam文件，应该是$wkd/clean/*bam~ 觉得自己超棒，有没有~
一天的时间，跑完了一个流程，但是，只跑了3个样本，服务器就over了，sad
明天再跟着洲更大神的视频跑一个流程，我是不是就算入门了呢？
（不交钱是我最后的倔强~）

附：jimmy原教程