看老大RNA-seq视频-P8/P9:RNA-seq:6-qc-

2019-01-01  本文已影响22人  小梦游仙境

P8/9:RNA-seq:6/7-qc质控

接下来要从fa转fq文件

ls -lh raw/|wc #看数据
less raw/nohup.out #有后台日志
grep failed raw/nohup.out
img
zless SRR1039508_1.fastq.gz #查看

非常长,要用报告来看,如何做报告,如下:

fastqc SRR1039510_1.fastq.gz
fastqc -t -10 SRR1039510_1.fastq.gz#-t -10加参数使速度加快
image-20181230123450921

(写不进去,是因为是老大的命令,要放进自己的路径)

ls *gz |xargs fastqc -t 10#批量生成fastqc报告

fastqc文件生产后,用multiqc,综合生产一个报告

multiqc ./
image-20181230165944176 image-20181230170009907 image-20181230170031425 image-20181230170048461 image-20181230170100528

老大视频里的未修改前的循环

image-20181230170510727

更改后的

image-20181230170554384

在后台运行

image-20181230170648457 image-20181230170827068

注:如果有config,在cat后改为config,如果没有就是$1

image-20181230172700029
source activate rna
bin_trim_galore=trim_galore
dir='/mm/project/airway/clean'
cat $1 |while read id
do
      arr=($id)
      fq1=${arr[0]}
      fq2=${arr[1]}
nohup $bin_trim_galore -q 25 --phred33 --length 36 -e 0.1 --stringency 3 --paired -o $dir $fq1 $fq2 &
done
source deactivate
dir='/mm/project/airway/clean'
cat config |while read id
do
      arr=($id)
      fq1=${arr[0]}
      fq2=${arr[1]}
nohup trim_galore -q 25 --phred33 --length 36 -e 0.1 --stringency 3 --paired -o $dir $fq1 $fq2 &
done

trim_galore一个

image-20181230182216099

之前就没跟着一起跑过这个循环,自己弄也跑不出来,和老大的代码一样也不行,改下也不行,真不知道怎么回事了.....就是没有我的任务😢

image-20181230185346564

config也没啥问题

image-20181230185530418

划重点了,我要默念一百遍!以上是我未跑出来的过程,为了记录过程,还是保留。最后还是我 老大一语中的,我明明是就改了dir的输出路径,出了问题,那我就应该把问题集中在路径这里就好了!

###这个是最后能跑的

source activate rna
bin_trim_galore=trim_galore
dir='/four/mm/project/airway/clean'  ###问题就在这个路径这,我没从根目录出发,是从mm这个相对路径出发的
cat $1 |while read id
do
      arr=($id)
      fq1=${arr[0]}
      fq2=${arr[1]}
nohup $bin_trim_galore -q 25 --phred33 --length 36 -e 0.1 --stringency 3 --paired -o $dir $fq1 $fq2 &
done
source deactivate
image-20181230221300389

已经跑了

image-20181230221405486

好像又出来什么问题,任务显示能跑,但是为什么是gzip

image-20181231001237200

单个跑

trim_galore --phred33 -q 25 -e 0.1 --length 36 --stringency 3 --paired -o ./ ~/fqmm/SRR1039508_1.fastq.gz ~/fqmm/SRR1039508_2.fastq.gz

P10:7-alignment-1

从老大视频里扒来的链接,先收藏

https://blog.csdn.net/xubo245/article/details/50878760

https://blog.csdn.net/xubo245/article/details/50836185

比对:hisat2、subjunc、star大多是针对转录组的,bwa、bowtie2是基因组

老大的参考基因组位置

/public/reference/genome/*

老大的参考基因组索引位置

/Public/reference/index/*

加了✳️和不加是不一样的,加了✳️可以把文件夹内的内容同时显示出来,不加的话只会显示当前路径下的文件夹

hisat有单独文件夹

/public/reference/index/hisat/*

每个文件取前1000行

ls ./*gz |while read ijfmklmf;do(zcat $ijfmklmf |head -1000 > 一个文件名);done

因为加了管道符,所以要加()

image-20181226153113419

由于输出到clean文件里了,所以想只要文件名

image-20181230200625914 image-20181230200654422 image-20181230201459975 image-20181230201619833

改名字

image-20181230202324986 image-20181230202702154

sam文件

image-20181230234856568 image-20181231143822782

我跟着做的

image-20181231153845434 image-20181231154315921 image-20181231154346749 image-20181231154654895

是不是链接过来的文件不能做下面这样的操作呢?

image-20181231160628179 image-20181231170134092 image-20181231170148669

上面是链接过来的文件就不行,我自己trim_galore的文件就可以,不知道为什么啊

image-20181231170526467

上面是链接过来的文件就不行,我自己trim_galore的文件就可以

image-20181231171723873

但是改成去掉“gz”的以后,就不能比对了😢哇,好神奇,知道了呢!是因为我按照老大的去掉了“.gz”后,就把fq.gz的文件给移动到别的文件夹了,然后就不可以了,因为RR1039508_1_val_1.fq只是SRR1039508_1_val_1.fq.gz的类似快捷方式的东西吧,并不是文件本身,-lh也能看到它没有大小。我在把fq.gz文件移动回来以后,就可以啦

image-20181231175921784

呃,但是报错了,是因为fq-2是0,这是为什么呢/😢

image-20181231180551941

视频里老大的比对:

image-20181231172848430

单个文件比对:

hisat2 -p 4 -x /four/mm/index/hisat/hg38/genome -1 SRR1039508_1_val_1.fq.gz -2 SRR1039508_2_val_2.fq.gz -S SRR1039508.tmp.sam#也可以输出为SRR1039508.hisat.sam,是没去掉gz的,要去掉gz的原因是有些软件会识别成压缩软件 

总之,不去掉.gz是可以比对的,就行了😄

用循环比对,hisat2比对到索引文件,然后samtools直接对生成的bam文件排序,以利于后续软件分析。

nohup cat SRR_Acc_List.txt |while read id;do  #复制一份id到当前路径下
hisat2 -p 5 -x ~/index/grch38/genome -1 \
${id}_1_val_1.fq.gz -2 ${id}_2_val_2.fq.gz | \
samtools sort -@ 5 -o ~/rna.GSE52778/sort.bam/${id}.sort.bam -
done &
image-20181231233220702 image-20181231231225867 image-20181231231330202

差异分析

崔老师ppt里的,也没跑出来

for fn in {508..523}
do
featureCounts -T 5 -p -t exon -g gene_id \
-a /four/mm/project/gtf/gencode.v29.annotation.gtf \ 
-o SRR1039$fn.counts.txt SRR1039$fn.sort.bam
# donot set dir
done

###我改完以后的,跑不出来
cat SRR_Acc_List.txt |while read id;do featureCounts -T 5 -p -t exon -g gene_id 
-a /four/mm/project/gtf/gencode.v29.annotation.gtf 
-o ${id}.counts.txt ${id}.sort.bam 
# donot set dir 
done
上一篇下一篇

猜你喜欢

热点阅读