脚本2：
perl /mnt/sdb/share/bin/deal_fa.pl
数据路径：/mnt/sdb/share/bin/linux-demo/program/fasta
要求：

1. 过滤掉长度1000bp一下的序列
2. 提取目标区域（region.txt）内的序列
3. 计算每个scaffold序列的长度
4. 以100bp为窗口计算scaffold窗口内的GC含量
5. scaffold序列基本信息统计（N50，N90，最大长度，最小长度等）
6. 以300个N链接scaffold序列
7. scaffold的拆分为contig
8. 每两个序列拆分为一个文件
9. 反向互补scaffold的序列
10. 把不标准的scaffold的序列变成标准的scaffold的序列
11. fasta文件转化为fastq格式文件
12. 修改fasta序列前缀名称；前缀为： part

【以下回答由 GPT 生成】

对于上述任务，你可以使用Perl编写一个脚本来完成。下面是一个可以实现这些功能的Perl脚本的示例代码：

```perl

!/usr/bin/perl

use strict; use warnings;

Task 1: 过滤掉长度低于1000bp的序列。

sub filter_sequences { my ($input_file, $output_file) = @_; open(my $in_fh, "<", $input_file) or die "Cannot open file: $!"; open(my $out_fh, ">", $output_file) or die "Cannot open file: $!"; my $sequence = ""; my $header = ""; while (my $line = <$in_fh>) { chomp($line); if ($line =~ /^>/) { if ($sequence ne "" && length($sequence) >= 1000) { print $out_fh "$header\n$sequence\n"; } $header = $line; $sequence = ""; } else { $sequence .= $line; } } if ($sequence ne "" && length($sequence) >= 1000) { print $out_fh "$header\n$sequence\n"; } close($in_fh); close($out_fh); }

Task 2: 从名为"region.txt"的文件中提取目标区域内的序列。

sub extract_sequences { my ($input_file, $region_file, $output_file) = @_; open(my $in_fh, "<", $input_file) or die "Cannot open file: $!"; open(my $region_fh, "<", $region_file) or die "Cannot open file: $!"; open(my $out_fh, ">", $output_file) or die "Cannot open file: $!"; my %regions; while (my $region_line = <$region_fh>) { chomp($region_line); my ($region_start, $region_end) = split("\t", $region_line); $regions{$region_start} = $region_end; } my $current_seq = ""; my $current_header = ""; while (my $line = <$in_fh>) { chomp($line); if ($line =~ /^>/) { if ($current_seq ne "" && exists $regions{$current_header}) { my $end_position = $regions{$current_header} > length($current_seq) ? length($current_seq) : $regions{$current_header}; print $out_fh "$current_header\n", substr($current_seq, 0, $end_position), "\n"; } $current_header = $line; $current_seq = ""; } else { $current_seq .= $line; } } if ($current_seq ne "" && exists $regions{$current_header}) { my $end_position = $regions{$current_header} > length($current_seq) ? length($current_seq) : $regions{$current_header}; print $out_fh "$current_header\n", substr($current_seq, 0, $end_position), "\n"; } close($in_fh); close($region_fh); close($out_fh); }

Task 3: 计算每个scaffold序列的长度。

sub calculate_lengths { my ($input_file) = @_; open(my $in_fh, "<", $input_file) or die "Cannot open file: $!"; my %lengths; my $current_header = ""; while (my $line = <$in_fh>) { chomp($line); if ($line =~ /^>/) { $current_header = $line; $lengths{$current_header} = 0; } else { $lengths{$current_header} += length($line); } } close($in_fh); return \%lengths; }

Task 4: 使用100bp的窗口计算每个scaffold序列窗口内的GC含量。

sub calculate_gc_content { my ($input_file, $window_size) = @_; open(my $in_fh, "<", $input_file) or die "Cannot open file: $!"; my %gc_content; my $current_header = ""; my $current_seq = ""; while (my $line = <$in_fh>) { chomp($line); if ($line =~ /^>/) { if ($current_header ne "") { for (my $i = 0; $i < length($current_seq); $i += $window_size) { my $window_sequence = substr($current_seq, $i, $window_size); my $gc_count = ($window_sequence =~ tr/GC//); my $window_start = $i + 1; my $window_end = $i + $window_size; push(@{$gc_content{$current_header}}, "$window_start-$window_end: " . ($gc_count / $window_size) * 100 . "%"); } } $current_header = $line; $current_seq = ""; } else { $current_seq .= $line; } } if ($current_header ne "") { for (my $i = 0; $i < length($current_seq); $i += $window_size) { my $window_sequence = substr($current_seq, $i, $window_size); my $gc_count = ($window_sequence =~ tr/GC//); my $window_start = $i + 1; my $window_end = $i + $window_size; push(@{$gc_content{$current_header}}, "$window_start-$window_end: " . ($gc_count / $window_size) * 100 . "%"); } } close($in_fh); return \%gc_content; }

Task 5: 统计scaffold序列的基本信息，包括N50、N90、最大长度和最小长度。

sub calculate_basic_stats { my ($input_file) = @_; open(my $in_fh, "<", $input_file) or die "Cannot open file: $!"; my @sequence_lengths; my $total_length = 0; my $current_length = 0; my $current_header = ""; my $max_length = 0; my $min_length; while (my $line = <$in_fh>) { chomp($line); if ($line =~ /^>/) { if ($current_length > 0) { push(@sequence_lengths, $current_length); $total_length += $current_length; if ($current_length > $max_length) { $max_length = $current_length; } if (!defined $min_length || $current_length < $min_length) { $min_length = $current_length; } } $current_header = $line; $current_length = 0; } else { $current_length += length($line); } } if ($current_length > 0) { push(@sequence_lengths, $current_length); $total_length += $current_length; if ($current_length > $max_length) { $max_length = $current_length; } if (!defined $min_length || $current_length < $min_length) { $min_length = $current_length; } } close($in_fh);

@sequence_lengths = sort { $b <=> $a } @sequence_lengths;
my $n50 = calculate_n(50, \@sequence_lengths, $total_length);
my $n90 = calculate_n(90, \@sequence_lengths, $total_length);

return ($max_length, $min_length, $n50, $n90);

}

Helper function for calculating N50 and N90

sub calculate_n { my ($percentage, $lengths, $total_length) = @_; my $threshold = $total_length * ($percentage / 100); my $current_length = 0; foreach my $length (@{$lengths}) { $current_length += $length; if ($current_length >= $threshold) { return $length; } } }

Task 6: 使用300个N将scaffold序列连接起来。

sub join_sequences { my ($input_file, $output_file) = @_; open(my $in_fh, "<", $input_file) or die "Cannot open file: $!"; open(my $out_fh, ">", $output_file) or die "Cannot open file: $!"; my $current_seq = ""; while (my $line = <$in_fh>) { chomp($line); if ($line =~ /^>/) { if ($current_seq ne "") { print $out_fh join("", ("N") x 300), "\n$current_seq\n"; } $current_seq = ""; } else { $current_seq .= $line; } } if ($current_seq ne "") { print $out_fh join("", ("N") x 300), "\n$current_seq\n"; } close($in_fh); close($out_fh); }

Task 7: 将scaffold序列拆分为contig。

sub split_sequences { my ($input_file) = @; open(my $in_fh, "<", $input_file) or die "Cannot open file: $!"; my $output_count = 1; my $output_file = "contig$output_count.fasta"; open(my $out_fh, ">", $output_file) or die "Cannot open file: $!"; while (my $line = <$in_fh>) { if ($line =~ /^>/) { close($out_fh); $output_count++; $output_file = "contig_$output_count.fasta"; open($out_fh, ">", $output_file) or die "Cannot open file: $!"; } print $out_fh $line; } close($in_fh); close($out_fh); }

Task 8: 将每两个序列拆分为一个文件。

sub split_files { my ($input_file) = @; open(my $in_fh, "<", $input_file) or die "Cannot open file: $!"; my $output_count = 1; my $current_header = ""; my $output_file = "output$output_count.fasta"; open(my $out_fh, ">", $output_file) or die "Cannot open file: $!"; while (my $line = <$in_fh>) { if ($line =~ /^>/) { if ($current_header ne "") { close($out_fh); $output_count++; $output_file = "output_$output_count.fasta"; open($out_fh, ">", $output_file) or die "Cannot open file: $!"; } $current_header = $line; } print $out_fh $line; } close($in_fh); close($out_fh); }

Task 9: 对scaffold序列进行反向互补操作。

sub reverse_complement { my ($input_file, $output_file) = @_; open(my $in_fh, "<", $input_file) or die "Cannot open file: $!"; open(my $out_fh, ">", $output_file) or die "Cannot open file: $!"; while (my $line = <$in_fh>) { chomp($line); if ($line =~ /^>/) { print $out_fh "$line\n"; } else { $line = reverse($line); $line =~ tr/ACGTacgt/TGCAtgca/; print $out_fh "$line\n"; } } close($in_fh); close($out_fh); }

Task 10: 将不标准的scaffold序列转换为标准的scaffold序列。

sub convert_sequences { my ($input_file, $output_file) = @; open(my $in_fh, "<", $input_file) or die "Cannot open file: $!"; open(my $out_fh, ">", $output_file) or die "Cannot open file: $!"; while (my $line = <$in_fh>) { chomp($line); if ($line =~ /^>/) { $line =~ s//-/g; } print $out_fh "$line\n"; } close($in_fh); close($out_fh); }

Task 11: 将FASTA文件转换为FASTQ格式文件。

sub convert_to_fastq { my ($input_file, $output_file) = @_; open(my $in_fh, "<", $input_file) or die "Cannot open file: $!"; open(my $out_fh, ">", $output_file) or die "Cannot open file: $!"; my $sequence = ""; my $header = ""; while (my $line = <$in_fh>) { chomp($line); if ($line =~ /^>/) { if ($sequence ne "") { $header =~ s/^>//; print $out_fh '@', $header, "\n", $sequence, "\n+\n", ("I" x length($sequence)), "\n"; } $header = $line; $sequence = ""; } else { $sequence .= $line; } } if ($sequence ne "") { $header =~ s/^>//; print $out_fh '@', $header, "\n", $sequence, "\n+\n", ("I" x length($sequence)), "\n"; } close($in_fh); close($out_fh); }

Task 12: 修改FASTA序列的前缀名称，将前缀改为"part"。

sub modify_prefix { my ($input_file, $output_file) = @; open(my $in_fh, "<", $input_file) or die "Cannot open file: $!"; open(my $out_fh, ">", $output_file) or die "Cannot open file: $!"; while (my $line = <$in_fh>) { chomp($line); if ($line =~ /^>/) { $line =~ s/^>.+?/>part_/; } print $out_fh "$line\n"; } close($in_fh); close($out_fh); }

Run the functions

Task 1: 过滤掉长度低于1000bp的序列。

filter_sequences("input.fasta", "filtered.fasta");

Task 2: 从名为"region.txt"的文件中提取目标区域内的序列。

extract_sequences("input.fasta", "region.txt", "extracted.fasta");

Task 3: 计算每个scaffold序列的长度。

my $lengths = calculate_lengths("input.fasta"); foreach my $header (keys %{ $lengths }) { print "Length of $header: ", $lengths->{$header}, "bp\n"; }

Task 4: 使用100bp的窗口计算每个scaffold序列窗口内的GC含量。

my $gc_content = calculate_gc_content("input.fasta", 100); foreach my $header (keys %{ $gc_content }) { print "GC Content of $header:\n"; foreach my $window (@{ $gc_content->{$header} }) { print "$window\n"; } }

Task 5: 统计scaffold序列的基本信息，包括N50、N90、最大长度和最小长度。

my ($max_length, $min_length, $n50, $n90) = calculate_basic_stats("input.fasta"); print "Max Length: $max_length\n"; print "Min Length: $min_length\n"; print "N50: $n50\n"; print "N90: $n90\n";

Task 6: 使用300个N将scaffold序列连接起来。

join_sequences("input.fasta", "joined.fasta");

Task 7: 将scaffold序列拆分为contig。

split_sequences("input.fasta");

Task 8: 将每两个序列拆分为一个文件。

split_files("input.fasta");

Task 9: 对scaffold序列进行反向互补操作。

reverse_complement("input.fasta", "reverse_complemented.fasta");

Task 10: 将不标准的scaffold序列转换为标准的scaffold序列。

convert_sequences("input.fasta", "converted.fasta");

Task 11: 将FASTA文件转换为FASTQ格式文件。

convert_to_fastq("input.fasta", "output.fastq");

Task 12: 修改FASTA序列的前