热门标签
热门文章
- 1next 14 appRouter redux数据持久化_nextjs redux 持久化存储
- 2mac与windows服务器 访问和共享
- 3Android 15全面解读：性能飙升、隐私守护与智能生活新纪元_安卓15
- 4标题：怎样通过Dialogflow构建一个聊天机器人？React版。_dialogflow机器人
- 5连接Mongodb数据库的步骤以及注意事项_如何连接mongodb数据库
- 6小程序公告php实现,小程序两种滚动公告栏的实现方法
- 7Git仓库完整迁移全过程_gitee 将a仓库的克隆到b仓库
- 8FP6381AS5CTR原厂SOT23-5 1.2A同步降压IC DC-DC变频器
- 9STM32参考代码，编译时出现“cannot open source input file, no such file or directory"错误
- 10微信小程序用户隐私保护指引设置指南_mp后台-设置-基本设置-服务内容声明-用户隐私保护指引]中声明“剪切板”隐私收集
当前位置: article > 正文
数据分析：宏基因组分析-从Raw data到profile的简单流程
作者：知新_RL | 2024-04-30 16:30:57
赞
踩
数据分析：宏基因组分析-从Raw data到profile的简单流程
简介
该宏基因组流程主要有四步，分别是1.检查raw data；2.获得高质量reads；3.合并PE数据；4. reads map到参考数据库得到profile。
初步想法：
- 先分别撰写每一步的基础脚本，如过滤，mapping等过程的脚本，只针对单样本；与此同时，设计好输入文件的格式；
- 接着脚本内部每个样本生成每个步骤的脚本，如sample1.trim.sh sample1.map.sh
- 然后将每步的脚本放置一起形成该步骤的综合脚本，如 step1.trim.sh
- 最后将含有每样本的各步骤的脚本综合在一起，为Run.all.sh
文件结构：脚本和结果文件
../MetaGenomics_pipeline/ ├── bin # 脚本 │ ├── humann.pl │ ├── kneaddata.pl │ ├── merge.pl │ ├── metaphlan.pl │ └── qc.pl ├── main.pl # 主程序 ├── result # 结果 │ ├── 00.quality │ │ ├── fastqc │ │ └── multiqc │ ├── 01.kneaddata │ ├── 02.merge │ ├── 03.humann │ │ ├── genefamilies │ │ ├── log │ │ ├── metaphlan │ │ ├── pathabundance │ │ └── pathcoverage │ ├── 04.metaphlan │ ├── Run.s1.qc.sh │ ├── Run.s2.kneaddata.sh │ ├── Run.s3.merge.sh │ ├── Run.s4.humann.sh │ ├── Run.s5.metaphlan.sh │ └── script # 每个样本的每一步脚本 │ ├── 00.quality │ ├── 01.kneaddata │ ├── 02.merge │ ├── 03.humann │ └── 04.metaphlan ├── Run.all.sh ├── test.tsv ├── TruSeq2-PE.fa -> /data/share/anaconda3/share/trimmomatic/adapters/TruSeq2-PE.fa └── work.sh # 启动脚本
- 1
- 2
- 3
- 4
- 5
- 6
- 7
- 8
- 9
- 10
- 11
- 12
- 13
- 14
- 15
- 16
- 17
- 18
- 19
- 20
- 21
- 22
- 23
- 24
- 25
- 26
- 27
- 28
- 29
- 30
- 31
- 32
- 33
- 34
- 35
- 36
步骤
先准备输入数据
find /RawData/ -name "*fq.gz" | sort | perl -e 'print "SampleID\tLaneID\tPath\n"; while(<>){chomp; $fq=(split("\/", $_))[-1]; $sampleid=$fq; $laneid=$fq; $sampleid=~s/\_R[1|2]\.fq.gz//g; $laneid=~s/\.fq.gz//g;print "$sampleid\t$laneid\t$_\n";}' > samples.fqpath.tsv
- 1
| SampleID | LaneID | Path |
|---|---|---|
| ND2 | ND2_R1 | RawData/ND2_R1.fq.gz |
| ND2 | ND2_R2 | RawData/ND2_R2.fq.gz |
| XL10 | XL10_R1 | RawData/XL10_R1.fq.gz |
| XL10 | XL10_R2 | RawData/XL10_R2.fq.gz |
| XL11 | XL11_R1 | RawData/XL11_R1.fq.gz |
| XL11 | XL11_R2 | RawData/XL11_R2.fq.gz |
| XL1 | XL1_R1 | RawData/XL1_R1.fq.gz |
| XL1 | XL1_R2 | RawData/XL1_R2.fq.gz |
| XL2 | XL2_R1 | RawData/XL2_R1.fq.gz |
Scan raw data
qc.pl： 使用fastqc和multiqc软件对raw data进行扫描，输入数据是 samples.fqpath.tsv，使用perl编程。
#!/usr/bin/perl use warnings; use strict; use Getopt::Long; my ($file, $real_dir, $out, $help, $version); GetOptions( "f|file:s" => \$file, "d|real_dir:s" => \$real_dir, "o|out:s" => \$out, "h|help:s" => \$help, "v|version" => \$version ); &usage if(!defined $out); # output my $dir_qc = "$real_dir/result/00.quality/fastqc"; system "mkdir -p $dir_qc" unless(-d $dir_qc); # script my $dir_script = "$real_dir/result/script/00.quality/"; system "mkdir -p $dir_script" unless(-d $dir_script); my @array_name; open(IN, $file) or die "can't open $file\n"; open(OT, "> $out") or die "can't open $out\n"; <IN>; while(<IN>){ chomp; my @tmp = split("\t", $_); if(-e $tmp[2]){ my $bash = join("", $dir_script, $tmp[1], ".fastqc.sh"); open(OT2, "> $bash") or die "can't open $bash\n"; print OT2 "fastqc -o $dir_qc --noextract $tmp[2]\n"; close(OT2); print OT "sh $bash\n"; } } close(IN); my $dir_mc = "$real_dir/result/00.quality/multiqc"; system "mkdir -p $dir_mc" unless(-d $dir_mc); print OT "multiqc $dir_qc --outdir $dir_mc\n"; close(OT); sub usage{ print <<USAGE; usage: perl $0 -f <file> -d <real_dir> -o <out> options: -f|file :[essential]. -d|real_dir :[essential]. -o|out :[essential]. USAGE exit; }; sub version { print <<VERSION; version: v1.1 update: 20201223 - 20201224 author: zouhua1\@outlook.com VERSION };
- 1
- 2
- 3
- 4
- 5
- 6
- 7
- 8
- 9
- 10
- 11
- 12
- 13
- 14
- 15
- 16
- 17
- 18
- 19
- 20
- 21
- 22
- 23
- 24
- 25
- 26
- 27
- 28
- 29
- 30
- 31
- 32
- 33
- 34
- 35
- 36
- 37
- 38
- 39
- 40
- 41
- 42
- 43
- 44
- 45
- 46
- 47
- 48
- 49
- 50
- 51
- 52
- 53
- 54
- 55
- 56
- 57
- 58
- 59
- 60
- 61
- 62
- 63
- 64
- 65
- 66
- 67
trim low quality reads and remove host DNA
kneaddata.pl：kneaddata内部自带过滤和比对软件
#!/usr/bin/perl use warnings; use strict; use Getopt::Long; my ($file, $real_dir, $out, $adapter, $help, $version); GetOptions( "f|file:s" => \$file, "d|real_dir:s" => \$real_dir, "a|adapt:s" => \$adapter, "o|out:s" => \$out, "h|help:s" => \$help, "v|version" => \$version ); &usage if(!defined $out); my $dir = "$real_dir/result/01.kneaddata"; system "mkdir -p $dir" unless(-d $dir); # script my $dir_script = "$real_dir/result/script/01.kneaddata/"; system "mkdir -p $dir_script" unless(-d $dir_script); open(IN, $file) or die "can't open $file"; my %file_name; <IN>; while(<IN>){ chomp; my @tmp = split("\t", $_); push (@{$file_name{$tmp[0]}}, $tmp[2]); } close(IN); my ($fq1, $fq2); open(OT, "> $out") or die "can't open $out\n"; foreach my $key (keys %file_name){ if (${$file_name{$key}}[0] =~ /R1/){ $fq1 = ${$file_name{$key}}[0]; }else{ $fq2 = ${$file_name{$key}}[0]; } if (${$file_name{$key}}[1] =~ /R2/){ $fq2 = ${$file_name{$key}}[1]; }else{ $fq1 = ${$file_name{$key}}[1]; } my $trim_opt = "ILLUMINACLIP:$adapter:2:40:15 SLIDINGWINDOW:4:20 MINLEN:50"; #if(-e $fq1 && -e $fq2){ my $bash = join("", $dir_script, $key, ".kneaddata.sh"); open(OT2, "> $bash") or die "can't open $bash\n"; print OT2 "kneaddata -i $fq1 -i $fq2 --output-prefix $key -o $dir -v -t 5 --remove-intermediate-output --trimmomatic /data/share/anaconda3/share/trimmomatic/ --trimmomatic-options \'$trim_opt\' --bowtie2-options \'--very-sensitive --dovetail\' -db /data/share/database/kneaddata_database/Homo_sapiens_Bowtie2_v0.1/Homo_sapiens\n"; close(OT2); print OT "sh $bash\n"; #} } print OT "kneaddata_read_count_table --input $dir --output $dir/01kneaddata_sum.tsv\n"; close(OT); sub usage{ print <<USAGE; usage: perl $0 -f <file> -d <real_dir> -o <out> -a <adapter> options: -f|file :[essential]. -d|real_dir :[essential]. -o|out :[essential]. -a|adapt :[essential]. USAGE exit; }; sub version { print <<VERSION; version: v1.1 update: 20201223 - 20201224 author: zouhua1\@outlook.com VERSION };
- 1
- 2
- 3
- 4
- 5
- 6
- 7
- 8
- 9
- 10
- 11
- 12
- 13
- 14
- 15
- 16
- 17
- 18
- 19
- 20
- 21
- 22
- 23
- 24
- 25
- 26
- 27
- 28
- 29
- 30
- 31
- 32
- 33
- 34
- 35
- 36
- 37
- 38
- 39
- 40
- 41
- 42
- 43
- 44
- 45
- 46
- 47
- 48
- 49
- 50
- 51
- 52
- 53
- 54
- 55
- 56
- 57
- 58
- 59
- 60
- 61
- 62
- 63
- 64
- 65
- 66
- 67
- 68
- 69
- 70
- 71
- 72
- 73
- 74
- 75
- 76
- 77
- 78
- 79
- 80
- 81
- 82
merge PE data
merge.pl：合并PE数据
#!/usr/bin/perl use warnings; use strict; use Getopt::Long; my ($file, $real_dir, $out, $help, $version); GetOptions( "f|file:s" => \$file, "d|real_dir:s" => \$real_dir, "o|out:s" => \$out, "h|help:s" => \$help, "v|version" => \$version ); &usage if(!defined $out); my $dir = "$real_dir/result/02.merge"; system "mkdir -p $dir" unless(-d $dir); # script my $dir_script = "$real_dir/result/script/02.merge/"; system "mkdir -p $dir_script" unless(-d $dir_script); open(IN, $file) or die "can't open $file"; my %file_name; <IN>; while(<IN>){ chomp; my @tmp = split("\t", $_); push (@{$file_name{$tmp[0]}}, $tmp[2]); } close(IN); my ($fq1, $fq2); open(OT, "> $out") or die "can't open $out\n"; foreach my $key (keys %file_name){ $fq1 = join("", "./result/01.kneaddata/", $key, "_paired_1.fastq"); $fq2 = join("", "./result/01.kneaddata/", $key, "_paired_2.fastq"); #if(-e $fq1 && -e $fq2){ my $bash = join("", $dir_script, $key, ".merge.sh"); open(OT2, "> $bash") or die "can't open $bash\n"; print OT2 "fastp -i $fq1 -I $fq2 -h $dir/$key\_merge.html -j $dir/$key\_merge.json -m --merged_out $dir/$key\_merge.fastq.gz --failed_out $dir/$key\_failed.fastq.gz --include_unmerged --overlap_len_require 6 --overlap_diff_percent_limit 20 --detect_adapter_for_pe -5 -r -l 20 -y --thread 5\n"; close(OT2); print OT "sh $bash\n"; #} } close(OT); sub usage{ print <<USAGE; usage: perl $0 -f <file> -d <real_dir> -o <out> options: -f|file :[essential]. -d|real_dir :[essential]. -o|out :[essential]. USAGE exit; }; sub version { print <<VERSION; version: v1.1 update: 20201223 - 20201224 author: zouhua1\@outlook.com VERSION };
- 1
- 2
- 3
- 4
- 5
- 6
- 7
- 8
- 9
- 10
- 11
- 12
- 13
- 14
- 15
- 16
- 17
- 18
- 19
- 20
- 21
- 22
- 23
- 24
- 25
- 26
- 27
- 28
- 29
- 30
- 31
- 32
- 33
- 34
- 35
- 36
- 37
- 38
- 39
- 40
- 41
- 42
- 43
- 44
- 45
- 46
- 47
- 48
- 49
- 50
- 51
- 52
- 53
- 54
- 55
- 56
- 57
- 58
- 59
- 60
- 61
- 62
- 63
- 64
- 65
- 66
- 67
- 68
- 69
get profile matrix
humann.pl：获取功能等profile数据
#!/usr/bin/perl use warnings; use strict; use Getopt::Long; my ($file, $real_dir, $out, $help, $version); GetOptions( "f|file:s" => \$file, "d|real_dir:s" => \$real_dir, "o|out:s" => \$out, "h|help:s" => \$help, "v|version" => \$version ); &usage if(!defined $out); my $dir = "$real_dir/result/03.humann"; system "mkdir -p $dir" unless(-d $dir); my $dir_log = "$real_dir/result/03.humann/log"; system "mkdir -p $dir_log" unless(-d $dir_log); my $genefamilies = "$real_dir/result/03.humann/genefamilies"; system "mkdir -p $genefamilies" unless(-d $genefamilies); my $pathabundance = "$real_dir/result/03.humann/pathabundance"; system "mkdir -p $pathabundance" unless(-d $pathabundance); my $pathcoverage = "$real_dir/result/03.humann/pathcoverage"; system "mkdir -p $pathcoverage" unless(-d $pathcoverage); my $dir_metaphlan = "$real_dir/result/03.humann/metaphlan"; system "mkdir -p $dir_metaphlan" unless(-d $dir_metaphlan); # script my $dir_script = "$real_dir/result/script/03.humann/"; system "mkdir -p $dir_script" unless(-d $dir_script); open(IN, $file) or die "can't open $file"; my %file_name; <IN>; while(<IN>){ chomp; my @tmp = split("\t", $_); push (@{$file_name{$tmp[0]}}, $tmp[2]); } close(IN); my ($fq); open(OT, "> $out") or die "can't open $out\n"; foreach my $key (keys %file_name){ $fq = join("", "./result/02.merge/", $key, "_merge.fastq.gz"); #if($fq){ my $bash = join("", $dir_script, $key, ".humann.sh"); open(OT2, "> $bash") or die "can't open $bash\n"; print OT2 "humann --input $fq --output $dir --threads 10\n"; print OT2 "mv $dir/$key\_merge_humann_temp/$key\_merge_metaphlan_bugs_list.tsv $dir_metaphlan/$key\_metaphlan.tsv\n"; print OT2 "mv $dir/$key\_merge_humann_temp/$key\_merge.log $dir_log\n"; print OT2 "mv $dir/$key\_merge_genefamilies.tsv $genefamilies/$key\_genefamilies.tsv\n"; print OT2 "mv $dir/$key\_merge_pathabundance.tsv $pathabundance/$key\_pathabundance.tsv\n"; print OT2 "mv $dir/$key\_merge_pathcoverage.tsv $pathcoverage/$key\_pathcoverage.tsv\n"; print OT2 "rm -r $dir/$key\_merge_humann_temp/\n"; close(OT2); print OT "sh $bash\n"; #} } close(OT); sub usage{ print <<USAGE; usage: perl $0 -f <file> -d <real_dir> -o <out> options: -f|file :[essential]. -d|real_dir :[essential]. -o|out :[essential]. USAGE exit; }; sub version { print <<VERSION; version: v1.1 update: 20201223 - 20201225 author: zouhua1\@outlook.com VERSION };
- 1
- 2
- 3
- 4
- 5
- 6
- 7
- 8
- 9
- 10
- 11
- 12
- 13
- 14
- 15
- 16
- 17
- 18
- 19
- 20
- 21
- 22
- 23
- 24
- 25
- 26
- 27
- 28
- 29
- 30
- 31
- 32
- 33
- 34
- 35
- 36
- 37
- 38
- 39
- 40
- 41
- 42
- 43
- 44
- 45
- 46
- 47
- 48
- 49
- 50
- 51
- 52
- 53
- 54
- 55
- 56
- 57
- 58
- 59
- 60
- 61
- 62
- 63
- 64
- 65
- 66
- 67
- 68
- 69
- 70
- 71
- 72
- 73
- 74
- 75
- 76
- 77
- 78
- 79
- 80
- 81
- 82
- 83
- 84
- 85
- 86
- 87
metaphlan.pl：获取物种组成谱
#!/usr/bin/perl use warnings; use strict; use Getopt::Long; my ($file, $real_dir, $out, $help, $version); GetOptions( "f|file:s" => \$file, "d|real_dir:s" => \$real_dir, "o|out:s" => \$out, "h|help:s" => \$help, "v|version" => \$version ); &usage if(!defined $out); my $dir = "$real_dir/result/04.metaphlan"; system "mkdir -p $dir" unless(-d $dir); # script my $dir_script = "$real_dir/result/script/04.metaphlan/"; system "mkdir -p $dir_script" unless(-d $dir_script); open(IN, $file) or die "can't open $file"; my %file_name; <IN>; while(<IN>){ chomp; my @tmp = split("\t", $_); push (@{$file_name{$tmp[0]}}, $tmp[2]); } close(IN); my ($fq); open(OT, "> $out") or die "can't open $out\n"; foreach my $key (keys %file_name){ $fq = join("", "./result/02.merge/", $key, "_merge.fastq.gz"); #if($fq){ my $bash = join("", $dir_script, $key, ".metaphlan.sh"); open(OT2, "> $bash") or die "can't open $bash\n"; print OT2 "metaphlan $fq --bowtie2out $dir/$key\_metagenome.bowtie2.bz2 --nproc 10 --input_type fastq -o $dir/$key\_metagenome.tsv --unknown_estimation -t rel_ab_w_read_stats\n"; close(OT2); print OT "sh $bash\n"; #} } print OT "merge_metaphlan_tables.py $dir/*metagenome.tsv > $dir/merge_metaphlan.tsv\n"; print OT "Rscript /data/share/database/metaphlan_databases/calculate_unifrac.R $dir/merge_metaphlan.tsv /data/share/database/metaphlan_databases/mpa_v30_CHOCOPhlAn_201901_species_tree.nwk $dir/unifrac_merged_mpa3_profiles.tsv"; close(OT); sub usage{ print <<USAGE; usage: perl $0 -f <file> -d <real_dir> -o <out> options: -f|file :[essential]. -d|real_dir :[essential]. -o|out :[essential]. USAGE exit; }; sub version { print <<VERSION; version: v1.0 update: 20201225 - 20201225 author: zouhua1\@outlook.com VERSION };
- 1
- 2
- 3
- 4
- 5
- 6
- 7
- 8
- 9
- 10
- 11
- 12
- 13
- 14
- 15
- 16
- 17
- 18
- 19
- 20
- 21
- 22
- 23
- 24
- 25
- 26
- 27
- 28
- 29
- 30
- 31
- 32
- 33
- 34
- 35
- 36
- 37
- 38
- 39
- 40
- 41
- 42
- 43
- 44
- 45
- 46
- 47
- 48
- 49
- 50
- 51
- 52
- 53
- 54
- 55
- 56
- 57
- 58
- 59
- 60
- 61
- 62
- 63
- 64
- 65
- 66
- 67
- 68
- 69
- 70
- 71
主程序
main.pl：生成所有的准备文件
#!/usr/bin/perl use warnings; use strict; use Getopt::Long; use FindBin qw($RealBin); use Cwd 'abs_path'; my ($file, $adapter, $out, $help, $version); GetOptions( "f|file:s" => \$file, "o|out:s" => \$out, "a|adapter:s" => \$adapter, "h|help:s" => \$help, "v|version" => \$version ); &usage if(!defined $out); my $Bin = $RealBin; my $cwd = abs_path; # output my $dir = "$cwd/result/"; system "mkdir -p $dir" unless(-d $dir); ########## output ######################################### # bash script per step # combine all steps in one script my $qc = join("", $dir, "Run.s1.qc.sh"); my $kneaddata = join("", $dir, "Run.s2.kneaddata.sh"); my $merge = join("", $dir, "Run.s3.merge.sh"); my $humann = join("", $dir, "Run.s4.humann.sh"); my $metaphlan = join("", $dir, "Run.s5.metaphlan.sh"); # scripts in bin my $bin = "$Bin/bin"; my $s_qc = "$bin/qc.pl"; my $s_kneaddata = "$bin/kneaddata.pl"; my $s_merge = "$bin/merge.pl"; my $s_humann = "$bin/humann.pl"; my $s_metaphlan = "$bin/metaphlan.pl"; ########## Steps in metagenomics pipeline ################# ################################################## # step1 reads quality scan `perl $s_qc -f $file -d $cwd -o $qc`; ################################################## # step2 filter and trim low quality reads; # remove host sequence `perl $s_kneaddata -f $file -d $cwd -a $adapter -o $kneaddata`; ################################################## # step3 merge PE reads `perl $s_merge -f $file -d $cwd -o $merge`; ################################################## # step4 get function profile `perl $s_humann -f $file -d $cwd -o $humann`; ################################################## # step5 get taxonomy profile `perl $s_metaphlan -f $file -d $cwd -o $metaphlan`; open(OT, "> $out") or die "can't open $out\n"; print OT "sh $qc\nsh $kneaddata\nsh $merge\nsh $humann\nsh $metaphlan\n"; close(OT); sub usage{ print <<USAGE; usage: perl $0 -f <file> -o <out> -a <adapter> options: -f|file :[essential]. -o|out :[essential]. -a|adapter :[essential]. USAGE exit; }; sub version { print <<VERSION; version: v1.1 update: 20201224 - 20201225 author: zouhua1\@outlook.com VERSION };
- 1
- 2
- 3
- 4
- 5
- 6
- 7
- 8
- 9
- 10
- 11
- 12
- 13
- 14
- 15
- 16
- 17
- 18
- 19
- 20
- 21
- 22
- 23
- 24
- 25
- 26
- 27
- 28
- 29
- 30
- 31
- 32
- 33
- 34
- 35
- 36
- 37
- 38
- 39
- 40
- 41
- 42
- 43
- 44
- 45
- 46
- 47
- 48
- 49
- 50
- 51
- 52
- 53
- 54
- 55
- 56
- 57
- 58
- 59
- 60
- 61
- 62
- 63
- 64
- 65
- 66
- 67
- 68
- 69
- 70
- 71
- 72
- 73
- 74
- 75
- 76
- 77
- 78
- 79
- 80
- 81
- 82
- 83
- 84
- 85
- 86
- 87
- 88
- 89
- 90
运行主程序：准备samples.fqpath.tsv和adapter.fa文件即可生成所有文件
perl main.pl -f samples.fqpath.tsv -a TruSeq2-PE.fa -o Run.all.sh
- 1
Run.all.sh 文件包含如下命令
sh /data/user/zouhua/pipeline/MetaGenomics_v2/result/Run.s1.qc.sh
sh /data/user/zouhua/pipeline/MetaGenomics_v2/result/Run.s2.kneaddata.sh
sh /data/user/zouhua/pipeline/MetaGenomics_v2/result/Run.s3.merge.sh
sh /data/user/zouhua/pipeline/MetaGenomics_v2/result/Run.s4.humann.sh
sh /data/user/zouhua/pipeline/MetaGenomics_v2/result/Run.s5.metaphlan.sh
- 1
- 2
- 3
- 4
- 5
声明：本文内容由网友自发贡献，不代表【wpsshop博客】立场，版权归原作者所有，本站不承担相应法律责任。如您发现有侵权的内容，请联系我们。转载请注明出处：https://www.wpsshop.cn/w/知新_RL/article/detail/514172
推荐阅读
相关标签
