写在前面
写了一个脚本,整理了一些命令。方便课题组所有人,具体功能是:
从Fastq格式的数据直接到Counts数据
下面介绍用法,
输入数据
- 也可以直接从NCBI下载数据之后进行转换(如果已经有原始数据,那么可以跳过)
首先是准备一个SRA编号的列表
vim SRAlist.txt
内容如下
SRR3080054
SRR3080055
SRR3080055
SRR3080057
SRR3080058
SRR3080059
使用以下命令下载
for file in `cat SRAlist.txt`;do perl -e '$id=shift;$prefix=substr($id,0,3);$run=substr($id,0,6);system qq{aria2c -j 20 ftp://ftp-trace.ncbi.nih.gov/sra/sra-instant/reads/ByRun/sra/$prefix/$run/$id/$id.sra}' $file;done
使用以下命令,进行数据转换,那么我们会得到1.样本数据文件
for file in *.sra;do
fasterq-dump --split-3 $file
done
1.样本数据
ls -1 S*.fastq
SRR3080055.sra_1.fastq
SRR3080055.sra_2.fastq
SRR3080056.sra_1.fastq
SRR3080056.sra_2.fastq
SRR3080057.sra_1.fastq
SRR3080057.sra_2.fastq
SRR3080058.sra_1.fastq
SRR3080058.sra_2.fastq
SRR3080059.sra_1.fastq
SRR3080059.sra_2.fastq
整理成一个表格
ls S*.fastq|paste - - > sample.list.txt
cat sample.list.txt
表格内容如下
SRR3080055.sra_1.fastq SRR3080055.sra_2.fastq
SRR3080056.sra_1.fastq SRR3080056.sra_2.fastq
SRR3080057.sra_1.fastq SRR3080057.sra_2.fastq
SRR3080058.sra_1.fastq SRR3080058.sra_2.fastq
SRR3080059.sra_1.fastq SRR3080059.sra_2.fastq
2.参考基因组
A_subgenome.fa
3.基因结构注释信息
A_subgenome.gff3
运行流程
perl /tools/XiaLabRNAseqPipe/Read2Counts.pl A_subgenome.fa sample.list A_subgenome.gff3
然后就等着就可以了
对于数据太多,硬盘空间有限
一次只运行10个样品,然后清理所有bam 和 bai
paraNum=10
ls *.fq.gz|paste - - |split -l $paraNum --additional-suffix .cc
for file in *.cc;do perl /tools/XiaLabRNAseqPipe/Read2Counts.pl /home/XiaLab/LabGenome/Litchi_Final_Genome/Lchinesis_genome.Chr.fasta $file /home/XiaLab/LabGenome/Litchi_Final_Genome/Lchinesis_genome.Chr.gff3;/bin/rm -r *.bam *.bai;done