3 wes测序质量的控制

原视频6:测序质量的控制首先建立文件夹$ cd ~/project/wes/$ mkdir {raw,clean,align,mutation,qc}这部分包括fastqc和multiqc两个软件查看测序质量,以及使用trim_galore软件进行过滤低质量reads和去除接头。1 QC1.1 fastqc没有原视频中文件,我用了下载的三个文件做例子。乳腺癌的组织样本。所以原视频中命令我也用不上,但是还是列出来find /public/project/wes/raw_data -name *.gz|grep -v '\._'|xargs fastqc -t 10 -o ./想知道为什么要-v ._,去看原视频中的文件命名。我的文件没那么复杂,可以下面这样$ find *.gz|xargs fastqc -t 20-t 20 一次运行20个文件。如果你有很多很多文件,参考我这篇批量对多个测序文件进行fastqc.1.2 multiqc$ multiqc -n wes ./......[INFO ] multiqc : This is MultiQC v1.0.dev0[INFO ] multiqc : Template : default[INFO ] multiqc : Searching './'[INFO ] fastqc : Found 8 reports[INFO ] multiqc : Report : wes.html[INFO ] multiqc : Data : wes_data[INFO ] multiqc : MultiQC complete结果如下所示:

multiqc结果

Per Sequence Quality Scores

接头以上结果发现,质量可以但是需要去接头2 trim-galore 过滤低质量reads和去接头Trim Galore! is a wrapper script to automate quality and adapter trimming as well as quality control#安装conda install -c bioconda trim-galorels /path/to/your/directory/*_1.fastq.gz >1ls /path/to/your/directory/*_2.fastq.gz >2paste 1 2 > config也可以用ls|grep >打开qc.sh,写入以下内容source activate wesbin_trim_galore=trim_galoredir='/mnt/f/kelly/bioTree/server/wesproject/raw_data'cat config |while read iddoarr=(${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运行上面脚本bash qc.sh就可以了。可以看到后台正在运行kelly@DESKTOP-MRA1M1F:~$ ps -efUID PID PPID C STIME TTY TIME CMDroot 1 0 0 17:06 ? 00:00:00 /initkelly 2 1 0 17:06 tty1 00:00:01 -bashkelly 527 1 5 22:46 tty1 00:00:27 perl /home/kell/miniconda3/bin/trim_galore -q 25 --phred33 --length 36kelly 528 1 5 22:46 tty1 00:00:30 perl /home/kell/miniconda3/bin/trim_galore -q 25 --phred33 --length 36kelly 529 1 5 22:46 tty1 00:00:29 perl /home/kell/miniconda3/bin/trim_galore -q 25 --phred33 --length 36kelly 530 1 5 22:46 tty1 00:00:26 perl /home/kell/miniconda3/bin/trim_galore -q 25 --phred33 --length 36运行时间大概4h,生成以下几个文件├── [ 88] 1├── [ 88] 2├── [ 176] config├── [ 50K] nohup.out├── [2.2G] SRR7696207_1.fastq.gz├── [4.8K] SRR7696207_1.fastq.gz_trimming_report.txt├── [2.1G] SRR7696207_1_val_1.fq.gz├── [2.4G] SRR7696207_2.fastq.gz├── [5.0K] SRR7696207_2.fastq.gz_trimming_report.txt├── [2.3G] SRR7696207_2_val_2.fq.gz├── [1.9G] SRR8517853_1.fastq.gz├── [4.8K] SRR8517853_1.fastq.gz_trimming_report.txt├── [1.8G] SRR8517853_1_val_1.fq.gz├── [2.1G] SRR8517853_2.fastq.gz├── [5.0K] SRR8517853_2.fastq.gz_trimming_report.txt├── [2.0G] SRR8517853_2_val_2.fq.gz├── [2.3G] SRR8517854_1.fastq.gz├── [4.8K] SRR8517854_1.fastq.gz_trimming_report.txt├── [2.2G] SRR8517854_1_val_1.fq.gz├── [2.6G] SRR8517854_2.fastq.gz├── [5.1K] SRR8517854_2.fastq.gz_trimming_report.txt├── [2.4G] SRR8517854_2_val_2.fq.gz├── [4.1G] SRR8517856_1.fastq.gz├── [4.9K] SRR8517856_1.fastq.gz_trimming_report.txt├── [4.0G] SRR8517856_1_val_1.fq.gz├── [4.5G] SRR8517856_2.fastq.gz├── [5.1K] SRR8517856_2.fastq.gz_trimming_report.txt├── [4.2G] SRR8517856_2_val_2.fq.gz└── [ 305] trim可见,各小了大概0.1G。其中,txt中的信息如下SUMMARISING RUN PARAMETERS==========================Input filename: SRR7696207_1.fastq.gzTrimming mode: paired-endTrim Galore version: 0.6.2Cutadapt version: 1.18Number of cores used for trimming: 1Quality Phred score cutoff: 25Quality encoding type selected: ASCII+33Adapter sequence: 'AGATCGGAAGAGC' (Illumina TruSeq, Sanger iPCR; auto-detected)Maximum trimming error rate: 0.1 (default)Minimum required adapter overlap (stringency): 3 bpMinimum required sequence length for both reads before a sequence pair gets removed: 36 bpOutput file will be GZIP compressedThis is cutadapt 1.18 with Python 2.7.15Command line parameters: -j 1 -e 0.1 -q 25 -O 3 -a AGATCGGAAGAGC SRR7696207_1.fastq.gzProcessing reads on 1 core in single-end mode ...Finished in 1011.58 s (36 us/read; 1.65 M reads/minute).=== Summary ===Total reads processed: 27,850,979Reads with adapters: 9,452,510 (33.9%)Reads written (passing filters): 27,850,979 (100.0%)Total basepairs processed: 4,177,646,850 bpQuality-trimmed: 21,999,885 bp (0.5%)Total written (filtered): 3,930,084,164 bp (94.1%)=== Adapter 1 ===Sequence: AGATCGGAAGAGC; Type: regular 3'; Length: 13; Trimmed: 9452510 times.No. of allowed errors:0-9 bp: 0; 10-13 bp: 1Bases preceding removed adapters:A: 21.7%C: 25.2%G: 28.1%T: 24.9%none/other: 0.0%Overview of removed sequenceslength count expect max.err error counts3 619264 435171.5 0 6192644 291812 108792.9 0 2918125 236788 27198.2 0 2367886 226036 6799.6 0 2260367 217969 1699.9 0 2179698 214073 425.0 0 2140739 230269 106.2 0 229591 67810 216909 26.6 1 209213 769611 223835 6.6 1 214662 917312 219837 1.7 1 210062 977513 219106 0.4 1 209421 968514 223045 0.4 1 212324 1072115 218505 0.4 1 208196 1030916 224812 0.4 1 213584 1122817 228422 0.4 1 216425 1199718 214056 0.4 1 204368 968819 216385 0.4 1 206368 1001720 207262 0.4 1 198505 875721 220284 0.4 1 209515 1076922 207937 0.4 1 198723 921423 203136 0.4 1 194604 853224 208015 0.4 1 198522 949325 200567 0.4 1 191904 866326 201338 0.4 1 192675 8663......

(0)

相关推荐

  • featureCounts和DEXseq做基于外显子定量的可变剪切

    rMATS这款差异可变剪切分析软件的使用体验 用LeafCutter探索转录组数据的可变剪切 用Expedition来分析单细胞转录组数据的可变剪切 使用SGSeq探索可变剪切 用DEXSeq分析可变 ...

  • LncRNA鉴定上游分析

    前面我们介绍了一系列的LncRNA鉴定相关文献的案例精选: 4个发育时间点的总共12个鸡转录组测序样本的长非编码RNA的鉴定 59匹马的8个组织的长非编码RNA的鉴定 9个组织的37个样本的大豆的长非 ...

  • 明码标价之免疫组库

    前面我带领大家通过IMGT数据库认知免疫组库,而且也一起从IMGT数据库下载免疫组库相关fasta序列,免疫组库重要的研究对象就是分成BCR的IGH,IGK,IGL这3类,以及TCR的TRA,TRB, ...

  • 使用阿里云+Docker分析RNA-Seq与ChIP-Seq

    写这篇文章的原因,是上一期大数据与生物信息结合的介绍 没有自己的服务器如何学习生物数据分析(上篇) 和 没有自己的服务器如何学习生物数据分析(下篇)(点击可以直接阅读)在生信技能树公众号发布以后,很多 ...

  • 使用gunzip命令的t参数检测fastq的gz文件完整度

    前面我们发布了 明码标价之普通转录组上游分析,终于开始接单了,第一个项目介绍98个转录组测序数据的表达量获取,超级简单,就是耗费计算资源,500G的fastq数据文件,中间步骤加起来,起码耗费2个T的 ...

  • 去rRNA可以解决GC含量双峰的右峰

    前些天我在生信技能树提出来了一个转录组数据分析的疑难杂症:RNA-seq的fastq文件里面为什么有gc含量的双峰,就是fastq测序数据质量控制的时候发现了GC含量的双峰,然后我简单分析了那些高重复 ...

  • ngs组学数据分析上下游分析都可以基于R语言吗?

    前些日子我们<生信技能树>的工程师做了一个ATAC-seq的项目,给客户汇报结果的时候,照例提供了全套代码.不过这次是从fq文件开始,所以大量代码都是在Linux平台的命令行而已,虽然给了 ...

  • 转录组比对步骤中的脚本运行问题与echo进行debug方法

    学徒在完成转录组实战任务的时候,出现了一个小bug,自顾自的纠结了一晚上不过来寻求我的帮助.让我非常生气,因为确实是很简单的两个小技巧,主要是理解全局变量.环境变量. 就是hisat2进行批量比对代码 ...

  • 跟着生信技能树学习microRNA测序学习

    我都忘记了自己的NGS流程系列增加了一个miRNA测序数据分析: 免费视频课程<RNA-seq数据分析> 免费视频课程<WES数据分析> 免费视频课程<ChIP-seq数 ...

  • 再次说明md5检查文件完整度的重要性

    最近服务器又停电,发现几个星期前提交的项目失败了几个样本: P5_DCIS  P2_Norm  P4_DCIS  P2_DCIS P9_DCIS  P10_Norm P9_Norm 所以我就去检查 c ...

  • 使用bowtie2去除宿主序列

    在研究组织或者肠道微生物时,常常需要去除宿主的DNA序列,以防止宿主的序列干扰研究.去宿主序列的主要研究方法是通过将质控后的序列与宿主基因组进行比对,将比对上的序列进行去除.比对软件通常有bowtie ...

  • 明码标价之普通转录组上游分析

    最近有粉丝在我们<生信技能树>公众号后台付费求助,想重新分析一下他自己领域的一个文章的转录组测序数据.因为那个文章并没有提供表达量矩阵,不过万幸的是人家上传了原始测序数据. 因为他的课题是 ...

  • lncRNA组装流程的软件介绍之trim-galore

    咱们<生信技能树>的B站有一个lncRNA数据分析实战,缺乏配套笔记,所以我们安排了100个lncRNA组装案例文献分享,以及这个流程会用到的100个软件的实战笔记教程! 下面是100个l ...