RNA-seq常用命令(无参)

简介: 0.前期准备先在工作目录下创建以下几个目录:01.raw_data #用于存放原始数据02.fastq #用于存放fastq格式数据03.

0.前期准备

先在工作目录下创建以下几个目录:

01.raw_data        #用于存放原始数据
02.fastq           #用于存放fastq格式数据
03.fastqc          #用于存放QC结果
04.trinity_result  #用于存放trinity结果
数据处理的顺序:01.raw_data > 02.fastq > 03.fastqc > 04.trinity_result

1.下载SRA

cd01.raw_data目录下,执行以下命令(二选一):

nohup axel -q ftp://ftp.ncbi.nlm.nih.gov/sra/sra-instant/reads/ByRun/sra/SRR/SRR130/SRR1300434/SRR1300434.sra &

axel是一种多线程下载工具,比wget快;-q参数能让axel进入静默模式(选用)。

nohup wget -c ftp://ftp.ncbi.nlm.nih.gov/sra/sra-instant/reads/ByRun/sra/SRR/SRR130/SRR1300434/SRR1300434.sra &

用wget的时候记得加-c参数,支持断点续传,就会方便很多。

用脚本批量下载数据的代码参考附录里的内容。


2.fastq-dump

  fastq-dump --defline-seq '@$sn[_$rn]/$ri' --split-files file.sra --split-3 -O|--outdir

这个是trinity给的命令,如果后面要用trinity进行无参组装的就需要加上--defline-seq '@$sn[_$rn]/$ri'这串看不懂的“乱码”,就要用这个命令,否则trinity会报错说无法识别这个格式。

经过改良后的命令:

nohup fastq-dump --defline-seq '@$sn[_$rn]/$ri' --split-3 file.sra -o ../02.fastq &

默认输出在上一层目录的02.fastq这个文件夹里。

-o:重定向输出结果的位置。将结果输出到上一级目录下的02.fastq文件夹内。

可选用的参数:

--gzip, --bzip2: 以压缩文件的方式输出结果。有利于减少文件的占用空间。


3.质控

不管什么数据,质控一下还是很有必要的,多看fastqc结果,即是积累经验,也是及时发现数据可能出现的问题。首先cd../02.fastq,之后执行:

  fastqc -t 8 *.fq -o ../03.fastqc

-t: 使用的线程数。现在用的服务器是20 CPU 40线程,以后这个值可以改大一点(range from 1~40,但是一般取10,如果文件比较多可以取20。)

-o: 重定向输出结果的位置。放在03.fastqc里面。

filezilla登陆服务器找到对应文件夹将生成的html文件下载到本地进行查看。如无问题则可进行下一步Trinity的操作,如有问题则进行质控部分余下的操作。

因为篇幅原因,将单独写一个关于如何查看fastqc结果的教程。敬请期待

需要去除接头的情况:

  • cutadapt

http://mp.weixin.qq.com/s/jZ37itt48slYrMNSvRxjNg

  • trimmomatic

https://www.jianshu.com/p/7b5591673255

需要对双端测序数据进行过滤低质量reads操作的情况:

  • sickle

https://mp.weixin.qq.com/s/ALY2FJCO10x02A7cy9qHyA

当你想改善你的数据的问题但是你又新入门啥软件都不会的时候(推荐使用):

  • fastp

https://www.jianshu.com/p/fae5b3c11e74

  fastp -i name.fastq -o name.fastp.fastq  #单端
  fastp -i name_1.fastq -o name_1.fastp.fastq -I name_2.fastq -O name_2.fastp.fastq  #双端

4.Trinity

cd../fastq文件夹后,执行以下命令:

双端测序命令:

nohup Trinity --seqType fq\fa --left ?  --right ?  --CPU 20 --max_memory 50G -output ../trinity_result/?.trinity &
  • left和right不能随意设置,left必须是_1,right必须是 _2.

  • 一般type是fq

  • ?的位置用具体需要操作的文件名代入即可。

单端测序命令:

nohup Trinity --seqType fq\fa --single ?  --CPU 20 --max_memory 50G -output ../trinity_result/?.trinity &

ps:输出结果的名字里必须含有trinity这个关键词,否则无法正常运行。


5.TransDecoder

cd../trinity_result/filename.trinity文件夹下后,找到一个叫做Trinity.fasta的文件,执行:

#step 1: 提取最长的开放阅读框
  TransDecoder.LongOrfs -t ./Trinity.fasta
  
  #step 2: (可选),BlastP搜索和Pfam搜索
    #BlastP搜索:蛋白库搜索, Swissprot (快) or Uniref90 (慢 but more comprehensive)
  blastp -query transdecoder_dir/longest_orfs.pep -db uniprot_sprot.fasta -max_target_seqs 1 -outfmt 6 -evalue 1e-5 -num_threads 10 > blastp.outfmt6
    #Pfam搜索:肽或蛋白域预测,需要安装hmmer3和Pfam数据库
  hmmscan --cpu 8 --domtblout pfam.domtblout /path/to/Pfam-A.hmm transdecoder_dir/longest_orfs.pep
  ​
  #step 3: 将Blast和Pfam搜索结果整合到编码区域选择这个步骤在尝试中发现比较费时间,建议用nohup命令进行操作。
  nohup TransDecoder.Predict -t target_transcripts.fasta --retain_pfam_hits pfam.domtblout --retain_blastp_hits blastp.outfmt6 &
  ​
  nohup TransDecoder.Predict -t trinity.fasta &  #简化版

通过TransDecoder这一步将能获得cds、pep、gff3、bed文件。用filezilla将以下五个文件拿到本地:

  • Trinity.fasta

  • Trinity.fasta.transdecoder.cds

  • Trinity.fasta.transdecoder.pep

  • Trinity.fasta.transdecoder.gff

  • Trinity.fasta.transdecoder.bed

    将这几个文件以可分辨的名字进行命名(如改成SRR1300215.fasta等)。


6.出个统计报告

还是在本目录下,对Trnity.fasta文件进行操作:

  TrinityStats.pl trinity_out_dir/Trinity.fasta >> stat.txt

会生成一个统计报告,并将内容重定向到stat.txt这个文件中。这个文件也记得保存下来。

结果展示。

以SRR1300215为例:

  ################################
  ## Counts of transcripts, etc.
  ################################
  Total trinity 'genes':  1524
  Total trinity transcripts:  1793
  Percent GC: 54.43
  ​
  ########################################
  Stats based on ALL transcript contigs:
  ########################################
  ​
    Contig N10: 1008
    Contig N20: 757
    Contig N30: 597
    Contig N40: 487
    Contig N50: 397
  ​
    Median contig length: 291
    Average contig: 379.04
    Total assembled bases: 679625
  ​
  ​
  #####################################################
  ## Stats based on ONLY LONGEST ISOFORM per 'GENE':
  #####################################################
  ​
    Contig N10: 965
    Contig N20: 715
    Contig N30: 559
    Contig N40: 461
    Contig N50: 382
  ​
    Median contig length: 287
    Average contig: 369.46
    Total assembled bases: 563055</pre>

7.附录

批量下载sra文件的脚本模板(mac版本,仅供参考)

  # 在终端运行
  #用axel下载会快一些
  brew install axel
  for i in {56..62}
  do
  # 也可以用wget 下载
  axel ftp://ftp-trace.ncbi.nlm.nih.gov/sra/sra-instant/reads/ByStudy/sra/SRP/SRP075/SRP075747/SRR35899$i/SRR35899$i.sra
  done

批量处理fastq-dump的脚本模板:(从人家的教程里摘下来的,仅供参考)

  #!/usr/bin/env bash
  for i in {56..62}
  do
  fastq-dump --gzip --split-3 -O /Users/chengkai/Desktop/zhuanlu_files -A SRR35899${i}.sra
  done

8.参考来源

文章:

微信公众号:

  • 生信菜鸟团

  • 生信者言


9.版权信息

Author:陈俊浩 Hanschen

QQ:1020607557

相关文章
|
数据挖掘
Seurat 4.0 | 单细胞转录组数据整合(scRNA-seq integration)
Seurat 4.0 | 单细胞转录组数据整合(scRNA-seq integration)
1766 0
Seurat 4.0 | 单细胞转录组数据整合(scRNA-seq integration)
|
4月前
|
数据可视化
Signac 单细胞|ATAC-seq Call peak
Signac 单细胞|ATAC-seq Call peak
47 0
Signac 单细胞|ATAC-seq Call peak
|
6月前
|
数据可视化 Java 数据处理
单细胞|RNA-seq & ATAC-seq 联合分析
单细胞|RNA-seq & ATAC-seq 联合分析
73 3
|
5月前
|
Linux Shell 数据处理
Linux命令seq的深入解析与应用
`seq`命令在Linux中用于生成数字序列,适用于数据处理和脚本编写。它支持自定义起始值、步长和结束值,可生成整数或浮点数。通过选项如`-f`(格式化输出)、`-s`(设置分隔符)、`-w`(宽度对齐)和`-r`(逆序)调整输出。在实践中,`seq`常与for循环结合,用于测试数据、文件命名等。注意增量为零会导致无限循环,格式和宽度选项不能同时使用。善用`seq`能提升命令行效率。
|
7月前
|
机器学习/深度学习 TensorFlow 算法框架/工具
seq2seq:中英文翻译
seq2seq:中英文翻译
57 1
|
7月前
|
存储 编解码 数据可视化
单细胞分析|整合 scRNA-seq 和 scATAC-seq 数据
单细胞分析|整合 scRNA-seq 和 scATAC-seq 数据
105 0
|
7月前
|
存储 数据可视化 数据挖掘
scRNA-seq|Seurat 整合分析
scRNA-seq|Seurat 整合分析
159 0
|
7月前
|
Linux
Linux命令(96)之seq
Linux命令(96)之seq
75 1
|
数据采集 算法 数据可视化
Seurat的scATAC-seq分析流程
Seurat的scATAC-seq分析流程
|
机器学习/深度学习 移动开发 自然语言处理
经典Seq2Seq与注意力Seq2Seq模型结构详解
经典Seq2Seq与注意力Seq2Seq模型结构详解
275 0
经典Seq2Seq与注意力Seq2Seq模型结构详解
下一篇
DataWorks