diff --git a/0_raw/step.sh b/0_raw/step.sh new file mode 100644 index 0000000..79106b7 --- /dev/null +++ b/0_raw/step.sh @@ -0,0 +1,9 @@ +###从NCBI下载基因组(DNA)/转录组(RNA)测序数据到获得预测基因的处理步骤(以Nyctotherus_ovalis物种为例,见/home/Data_01/yizhenzhen06/public_share/Saccharomyces_cerevisiae_DNA。如果没有注明仅基因组数据要做,其他步骤基因组和转录组步骤一样。直接写在本文件里的美丽表示可以直接终端运行指令,其他指令要提交pbs任务运行。 + +screen -S ncbi_dl #screen命令创建一个新的窗口(这里窗口命名为ncbi_dl)用来下载原始测序数据 +wget -c -b -i species_name_urls.txt #NCBI下载每个物种相应所有原始数据至服务器,species_name_urls.txt包括SRA数据的urls链接 +#命令运行成功后,按下组合键ctrl+A+D退出当前窗口 +screen -r ncbi_dl #返回ncbi_dl窗口,查看下载进度 + +:fastq-dump --defline-seq '@$sn[_$rn]/$ri' --split-files SRR10355986.1 #RNA数据转换数据格式为左右两端read文件 +:fastq-dump --split-files SRR10355986.1 #DNA数据转换数据格式为左右两端read文件 diff --git a/1_trim/step.sh b/1_trim/step.sh new file mode 100644 index 0000000..5dd717d --- /dev/null +++ b/1_trim/step.sh @@ -0,0 +1,23 @@ +###从NCBI下载基因组(DNA)/转录组(RNA)测序数据到获得预测基因的处理步骤(以Nyctotherus_ovalis物种为例,见/home/Data_01/yizhenzhen06/public_share/Saccharomyces_cerevisiae_DNA。如果没有注明仅基因组数据要做,其他步骤基因组和转录组步骤一样。直接写在本文件里的美丽表示可以直接终端运行指令,其他指令要提交pbs任务运行。 + +fastqc -o 1_trim/QC/ --noextract -t 2 -f fastq 0_rawdata/SRRXXXXXX_1.fastq 0_rawdata/SRRXXXXXX_2.fastq + +#分别检查左、右端数据质量,数据量大可以一起提交脚本运行。 +#如果数据质量不好,就做trimmomatic质控过滤。如果质量没问题,就分别合并所有左、右端两端reads文件进行组装。 + +:cat SRRXXXXXX.1_1.fastq SRRXXXXXX.1_1.fastq > species_name_1.fastq #合并所有左端reads文件 +:cat SRRXXXXXX.1_2.fastq SRRXXXXXX.1_2.fastq > species_name_2.fastq #合并所有右端reads文件 + +trimmomatic PE -phred33 -threads 16 \ +1_trim/species_name_1.fastq \ +1_trim/species_name_2.fastq \ +1_trim/species_name_1_p.fastq \ +1_trim/species_name_1_up.fastq \ +1_trim/species_name_2_p.fastq \ +1_trim/species_name_2_up.fastq \ +ILLUMINACLIP:../../adapters/NexteraPE-PE.fa:2:30:10:8:true LEADING:3 TRAILING:3 SLIDINGWINDOW:5:20 MINLEN:36 HEADCROP:15 +#trimmomatic质控过滤,参数需根据具体数据质量进行调整,大型数据请提交脚本运行。 + +fastqc -o 1_trim/QC/ --noextract -t 2 -f fastq 0_rawdata/SRRXXXXXX_1.fastq 0_rawdata/SRRXXXXXX_2.fastq + +#再次检查左、右端数据质量,数据量大同样提交脚本运行。 \ No newline at end of file diff --git a/2_assembly/step.sh b/2_assembly/step.sh new file mode 100644 index 0000000..4515f18 --- /dev/null +++ b/2_assembly/step.sh @@ -0,0 +1,18 @@ +###从NCBI下载基因组(DNA)/转录组(RNA)测序数据到获得预测基因的处理步骤(以Nyctotherus_ovalis物种为例,见/home/Data_01/yizhenzhen06/public_share/Saccharomyces_cerevisiae_DNA。如果没有注明仅基因组数据要做,其他步骤基因组和转录组步骤一样。直接写在本文件里的美丽表示可以直接终端运行指令,其他指令要提交pbs任务运行。 + +Trinity --seqType fq --max_memory 50G --no_version_check \ +--left species_name_1_p.fastq \ +--right species_name_2_p.fastq \ +--CPU XX +#RNA数据用Trinity组装,--CPU调整调用cpu核心数,其它参数需根据具体数据情况进行调整,提交脚本运行。 + +spades.py --sc -t XX -m 230 \ +-1 species_name_1_p.fastq \ +-2 species_name_2_p.fastq \ +-o spades_species_name +#DNA数据用spades组装,-t调整调用cpu核心数,其它参数需根据具体数据情况进行调整,提交脚本运行。 + +busco -i spades/scaffolds.fasta -m XXX \ +-l /home/Data_01/share_databases/busco/alveolata_odb10/alveolata_odb10 \ +-o busco_out_alve --offline -f +#BUSCO 检测完整度,-m调整评估模式(trans转录组 geno基因组) \ No newline at end of file diff --git a/3_decontamination/step.sh b/3_decontamination/step.sh new file mode 100644 index 0000000..70111b4 --- /dev/null +++ b/3_decontamination/step.sh @@ -0,0 +1,58 @@ +###从NCBI下载基因组(DNA)/转录组(RNA)测序数据到获得预测基因的处理步骤(以Nyctotherus_ovalis物种为例,见/home/Data_01/yizhenzhen06/public_share/Saccharomyces_cerevisiae_DNA。如果没有注明仅基因组数据要做,其他步骤基因组和转录组步骤一样。直接写在本文件里的美丽表示可以直接终端运行指令,其他指令要提交pbs任务运行。 + +blastn -db /home/Data_01/share_databases/ncbi_genome_refseq_bac_old/ncbi_genome_refseq_bac_old \ +-query /home/Data_01/yizhenzhen06/public_share/Saccharomyces_cerevisiae_DNA/2_assembly/spades/scaffolds.fasta \ +-out 1_Sacc_cerevisiae_dec_bac.tab \ +-outfmt "6 qseqid qlen sseqid slen pident qcovs qstart qend sstart send mismatch gapopen evalue bitscore" \ +-evalue 1e-5 -max_target_seqs 1 -num_threads X + +sed -i '1i\tquery ID\tquery sequence length\tsubject ID\tsubject qequence length\tpercent identities\tcoverage\tqstart\tqend\tsstart\tsend\tmismatch\tgapopen\tevalue\tbitscores ' 1_Sacc_cerevisiae_dec_bac.tab +#去细菌基因组污染,两条指令放一个脚本运行(注意检查数据库路径) + +awk '{ if ($5 > 80 && $6 > 70) print $1 }' 1_Sacc_cerevisiae_dec_bac.tab > 1_Sacc_cerevisiae_dec_bac.id #提取第5列相似度大于80%且覆盖率大于70%的序列ID(这一步建议需要下载到本地手动检查输出的id条数是否和Sacc_cerevisiae_dec_bac.tab(手动excel筛选)文件的第5列相似度大于80%且覆盖率大于70%的序列条数一致) +seqkit grep -v -f Sacc_cerevisiae_dec_bac.id ../2_assembly/trinity_out_dir/Trinity.fasta > 1_Sacc_cerevisiae_dec_bac.fasta #从组装结果文件中删去细菌污染序列 + +#################################################################################################### + +blastn -db /home/Data_01/share_databases/human_genome_new/GCA_021234545.1_HCC1395BL_v1.0_genomic_db \ +-query 1_Sacc_cerevisiae_dec_bac.fasta \ +-out 2_Sacc_cerevisiae_dec_bac_hum.tab \ +-outfmt "6 qseqid qlen sseqid slen pident qcovs qstart qend sstart send mismatch gapopen evalue bitscore" \ +-evalue 1e-5 -max_target_seqs 1 -num_threads X + +sed -i '1i\tquery ID\tquery sequence length\tsubject ID\tsubject qequence length\tpercent identities\tcoverage\tqstart\tqend\tsstart\tsend\tmismatch\tgapopen\tevalue\tbitscores ' 2_Sacc_cerevisiae_dec_bac_hum.tab +#去人类基因组污染,两条指令放一个脚本运行(注意检查数据库路径) + +awk '{ if ($5 > 80 && $6 > 70) print $1 }' 2_Sacc_cerevisiae_dec_bac_hum.tab > 2_Sacc_cerevisiae_dec_bac_hum.id #提取第5列相似度大于80%且覆盖率大于70%的序列ID(这一步建议需要下载到本地手动检查输出的id条数是否和2_Sacc_cerevisiae_dec_bac_hum.tab(手动excel筛选)文件的第5列相似度大于80%且覆盖率大于70%的序列条数一致) +seqkit grep -v -f 2_Sacc_cerevisiae_dec_bac_hum.id 1_Sacc_cerevisiae_dec_bac.fasta > 2_Sacc_cerevisiae_dec_bac_hum.fasta #从组装结果文件中删去人类基因组污染序列 + +#################################################################################################### + +blastn -db /home/Data_01/share_databases/human_genome_new/GCF_000001405.39_GRCh38.p13_rna_db \ +-query 2_Sacc_cerevisiae_dec_bac_hum.fasta \ +-out 3_Sacc_cerevisiae_dec_bac_hum_geno_rna.tab \ +-outfmt "6 qseqid qlen sseqid slen pident qcovs qstart qend sstart send mismatch gapopen evalue bitscore" \ +-evalue 1e-5 -max_target_seqs 1 -num_threads X + +sed -i '1i\tquery ID\tquery sequence length\tsubject ID\tsubject qequence length\tpercent identities\tcoverage\tqstart\tqend\tsstart\tsend\tmismatch\tgapopen\tevalue\tbitscores ' 3_Sacc_cerevisiae_dec_bac_hum_geno_rna.tab +#去人类转录组污染,两条指令放一个脚本运行(注意检查数据库路径) + +awk '{ if ($5 > 80 && $6 > 70) print $1 }' 3_Sacc_cerevisiae_dec_bac_hum_geno_rna.tab > 3_Sacc_cerevisiae_dec_bac_hum_geno_rna.id #提取第5列相似度大于80%且覆盖率大于70%的序列ID(这一步需要下载到本地手动检查输出的id条数是否和Nyct_bac_human_genome_blastn.tab(excel筛选)文件的第5列相似度大于80%且覆盖率大于70%的序列条数一致) +seqkit grep -v -f 3_Sacc_cerevisiae_dec_bac_hum_geno_rna.id 2_Sacc_cerevisiae_dec_bac_hum.fasta > 3_Sacc_cerevisiae_dec_bac_hum_geno_rna.fasta #从组装结果文件中删去人类转录组污染序列 + +#################################################################################################### + +seqkit seq -n -i ../2_assembly/spades/scaffolds.fasta > 4_Sacc_cerevisiae_dec_all.txt +#提取序列ID,仅基因组数据要做 +awk '{if($1~"_0.")print$1}' 4_Sacc_cerevisiae_dec_all.txt > 4_Sacc_cerevisiae_dec_all_cov1.id +#提取coverage小于1的序列ID,仅基因组数据要做 +seqkit grep -v -f 4_Sacc_cerevisiae_dec_all_cov1.id 4_Sacc_cerevisiae_dec_all.fasta > 4_Sacc_cerevisiae_dec_all_cov1.fasta +#去除coverage小于1的序列,仅基因组数据要做 + +#################################################################################################### + +cd-hit-est \ +-i 4_Sacc_cerevisiae_dec_all_cov1.fasta \ +-o 5_Sacc_cerevisiae_dec_all_cov1_cdhit.fasta \ +-c 0.9 -n 8 -d 0 -M 32000 -T 48 +#去冗余(核酸序列) \ No newline at end of file diff --git a/4_prediction/step.sh b/4_prediction/step.sh new file mode 100644 index 0000000..4ef590a --- /dev/null +++ b/4_prediction/step.sh @@ -0,0 +1,34 @@ +###从NCBI下载基因组(DNA)/转录组(RNA)测序数据到获得预测基因的处理步骤(以Nyctotherus_ovalis物种为例,见/home/Data_01/yizhenzhen06/public_share/Saccharomyces_cerevisiae_DNA。如果没有注明仅基因组数据要做,其他步骤基因组和转录组步骤一样。直接写在本文件里的美丽表示可以直接终端运行指令,其他指令要提交pbs任务运行。 + +seqkit sample ../3_decontamination/5_Sacc_cerevisiae_dec_all_cov1_cdhit.fasta \ +-o 5_Sacc_cerevisiae_dec_all_cov1_cdhit_sample.fasta \ +-p 0.25 -s 11 +#seqkit抽提 1/2 或 1/4 序列 -p 参数调整 + +#Codetta 读取基因组,然后利用已知蛋白质的数据库,通过比对保守蛋白质序列的隐马尔可夫模型(Profile HMMs)与输入序列,计算出每个密码子最可能对应的氨基酸 +#如果一个基因编码某种蛋白质,该蛋白质的氨基酸序列应该与已知蛋白质具有一定的相似性,Codetta 就通过寻找这种相似性来推断每个密码子的氨基酸解码 + +codetta.py 5_Sacc_cerevisiae_dec_all_cov1_cdhit_sample.fasta +#codetta预测密码子表(进入 py36 环境,提交任务) + +#################################################################################################### + +augustus \ +--progress=true --protein=on --introns=on --start=on --stop=on --cds=on --codingseq=on --gff3=on \ +--species=tetrahymena \ +--extrinsicCfgFile=/home/Data_01/yizhenzhen06/miniconda3/envs/augustus/config/extrinsic/extrinsic.M.RM.E.W.cfg \ +../3_decontamination/5_Sacc_cerevisiae_dec_all_cov1_cdhit.fasta \ +--outfile=Sacc_cerevisiae.gff3 +#基因组数据基因预测,提交脚本运行 +getAnnoFasta.pl Sacc_cerevisiae.gff3 && mv Sacc_cerevisiae.aa Sacc_cerevisiae_augustus.fasta +# 从预测结果中提取氨基酸序列 + +TransDecoder.LongOrfs \ +-t ../3_decontamination/5_Sacc_cerevisiae_dec_all_cov1_cdhit.fasta \ +-O LongOrfs/Sacc_cerevisiae_clean_cdhit_ciliate \ +-G Ciliate +TransDecoder.Predict \ +-t ../3_decontamination/5_Sacc_cerevisiae_dec_all_cov1_cdhit.fasta \ +-O Sacc_cerevisiae_clean_cdhit_predict \ +-G Ciliate +#转录组数据基因预测,-G参数改成help命令中相应物种的对应密码子,两条指令一起提交脚本运行 \ No newline at end of file