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
从老大视频里扒来的链接，先收藏
比对：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
网友评论