【Linux】bio_funclib函数库
Dec 12, 2023
»
plogs
参考：
- gibbslab/biobash: Common BASH routines for computational biology (github.com)
- Linux练习册07.3-Shell脚本在生信中的应用-知识点.pdf
函数库总架构
个人认为一个sh文件做一件事情比较好…所以这里的每个“函数”都是一个独立的可执行文件，最后将这些文件整合在一起，做了一个bio_funclib函数库。
| 函数 | 注释 |
|---|---|
| statBioFile | 给一个文件（可能是fa, fq, BED, GFF或其它），对文件进行简单统计，输出统计结果 |
| fq2fa | 将fq转fa文件，用户提供fq格式的文件，将其转换为fa格式 |
| revCmpDNAseq | 对DNA序列进行反转互补，用户可以选择仅反转或仅互补，或反转且互补 |
| DNA2RNA | 将DNA序列转化为RNA序列，但也可以选择RNA转DNA |
0. 示例文件介绍
| 文件名 | 文件格式 |
|---|---|
| Homo_sapiens.GRCh38.110.abinitio.gff3 | gff3 |
| homo_sapiens.GRCh38.A549.Regulatory_Build.regulatory_activity.20221007.gff | gff |
| Homo_sapiens.GRCh38.cdna.abinitio.fa | fasta |
| sarsCov2Genomes.fasta | fasta |
| SRR3091420_1_chr6.fastq | fastq |
| chr.bed | bed |
1. statBioFile
给一个文件（可能是fa, fq, BED, GFF或其它），对文件进行简单统计，输出统计结果。
- 对于fa文件，输出
- seq_count - fasta序列长度
- base_count - 四种碱基数
- gc_count - GC含量%
- length_range - 序列的最长和最短长度
- 对于fq文件，输出
- base_total - 碱基总数
- n_reads - 含有 N 碱基的reads个数
- n_count - 所有的reads中N碱基的总数
- ATCGNdis - 所有序列的第一位碱基的ATCGNatcg分布情况
- 对于bed文件，输出
- lines - bed文件的行数
- 对于gff文件，输出
- columns - gff文件的行数
- 对于不在此函数中的文件，如txt文件，提示
"Unsupported file type: $fileType"
#!/bin/bash
# 帮助信息
usage() {
cat <<HELP
Usage: $(basename $0) <file>
The script will perform statistics on various biological data files.
Supported File Types:
.fa fasta file - outputs sequence count, base count, GC content, and length range.
.fq fastq file - outputs base total, number of reads with N base, total N base count, and distribution of the first base.
.bed BED file - outputs number of lines, and start and end positions.
.gff GFF file - outputs number of columns.
For unsupported file types, the script will provide an appropriate message.
HELP
exit 1
}
handle_file() {
local file_path=$1
local file_extension="${file_path##*.}"
local file_type=$(echo "$file_extension" | awk '{print tolower($0)}')
case $file_type in
fa | fasta)
echo "Processing fasta file: $file_path"
local seq_count=$(grep -c '^>' "$file_path")
local base_count=$(grep -v '^>' "$file_path" | grep -o . | wc -l)
local gc_count=$(grep -v '^>' "$file_path" | grep -o '[GC]' | wc -l)
local length_ranges=$(grep -v '^>' "$file_path" | awk '{print length}' | sort -n | awk 'NR==1; END{print}')
echo "* Sequence count: $seq_count"
echo "* Base count: $base_count"
echo "* GC count: $gc_count"
echo "* GC content: $(echo "scale=2; ($gc_count/$base_count)*100" | bc)%"
echo "* Length range: $length_ranges"
;;
fq | fastq )
echo "Processing fastq file: $file_path"
local base_total=$(awk 'NR%4==2' "$file_path" | wc -c)
local n_reads=$(awk 'NR%4==2' "$file_path" | grep -c 'N')
local n_count=$(awk 'NR%4==2' "$file_path" | grep -o 'N' | wc -l)
local atcgn_distribution=$(awk 'NR%4==2' "$file_path" | cut -c1 | sort | uniq -c)
echo "* Base total: $base_total"
echo "* Reads with N: $n_reads"
echo "* N count: $n_count"
echo "* ATCGN distribution: $atcgn_distribution"
;;
bed)
echo "Processing BED file: $file_path"
local lines=$(wc -l < "$file_path")
echo "* Number of lines: $lines"
;;
gff)
echo "Processing GFF file: $file_path"
local columns=$(awk '{print NF; exit}' "$file_path")
echo "* Number of columns: $columns"
;;
*)
echo "Unsupported file type: .$file_type"
return 1
;;
esac
return 0
}
# 检查是否存在至少一个参数
if [ "$#" -eq 0 ]; then
usage
fi
file_path="$1"
# 检查文件路径是否为空以及文件是否存在
if [ -z "$file_path" ]; then
echo "Error: No file path provided."
usage
elif [ ! -f "$file_path" ]; then
echo "Error: File does not exist."
usage
else
# Call the generic function to handle the file
handle_file "$file_path"
fi
使用说明：
# 赋予可执行权限
chmod +x statBioFile
# 使用脚本
./statBioFile /path/to/your/datafile.fa
运行示例：
➜ ./statBioFile Homo_sapiens.GRCh38.cdna.abinitio.fa
Processing fasta file: Homo_sapiens.GRCh38.cdna.abinitio.fa
* Sequence count: 50174
* Base count: 62676097
* GC count: 33561747
* GC content: 53.00%
* Length range: 1
60
➜ ./statBioFile sarsCov2Genomes.fasta
Processing fasta file: sarsCov2Genomes.fasta
* Sequence count: 17
* Base count: 508282
* GC count: 183308
* GC content: 36.00%
* Length range: 14
60
➜ ./statBioFile SRR3091420_1_chr6.fastq
Processing fastq file: SRR3091420_1_chr6.fastq
* Base total: 41758450
* Reads with N: 9975
* N count: 11654
* ATCGN distribution: 131514 A
377795 C
262324 G
519 N
63017 T
➜ ./statBioFile homo_sapiens.GRCh38.A549.Regulatory_Build.regulatory_activity.20221007.gff
Processing GFF file: homo_sapiens.GRCh38.A549.Regulatory_Build.regulatory_activity.20221007.gff
* Number of columns: 12
➜ ./statBioFile chr.bed
Processing BED file: chr.bed
* Number of lines: 16959
2. fq2fa
将fq转fa文件，用户提供fq格式的文件，将其转换为fa格式。
在FASTQ文件中，每个序列由四行表示：
- 以’@’开头的行，包含序列标识符和可选的描述（类似于FASTA文件的标题行）。
- 序列本身。
- 以’+’开头的行，可选地跟随相同的序列标识符和描述。
- 包含序列质量分数的行。
而在FASTA文件中，每个序列由两行表示：
- 以’>’开头的行，包含序列标识符和可选的描述。
- 序列本身。
当从FASTQ格式转换为FASTA格式时，实际上是在移除以’+’开头的行和质量分数行，这样总行数减少了一半。因此，转换后的FASTA文件的行数应该是原始FASTQ文件行数的一半。可以利用这一点来检查确认转换过程如预期进行。
#!/bin/bash
# 检查是否存在至少一个参数
if [ "$#" -ne 1 ]; then
echo "Usage: $0 <input_fastq_file>"
exit 1
fi
input_fastq_file="$1"
output_fasta_file="${input_fastq_file%.fq}.fa"
# 这里是用的Linux练习册中的命令行
awk '{if (NR%4==1||NR%4==2) print}' "$input_fastq_file" | sed 's/@/>/g' > "$output_fasta_file"
# 检查行数，以免出错
echo "Check of line counts:"
fq_line_count=$(cat "$input_fastq_file" | wc -l)
fa_line_count=$(sed -n '1~4p;2~4p' "$input_fastq_file" | wc -l)
echo "The original fastq file has $fq_line_count lines."
echo "The converted fasta file should have approximately half this number: $fa_line_count lines."
使用说明：
# 赋予可执行权限
chmod +x fq2fa
# 使用脚本
./fq2fa /path/to/your/datafile.fq
但是当我使用这个脚本时，有一个小的报错：
sed: 1: "1~4p;2~4p": invalid command code ~
原因是：~ 运算符在sed中用于指定在每n行上操作的步长值，但在非Linux系统 (如macOS) 上可能不兼容（我就是macOS系统）。更改后的代码如下。
#!/bin/bash
# 检查是否存在至少一个参数
if [ "$#" -ne 1 ]; then
echo "Usage: $0 <input_fastq_file>"
exit 1
fi
input_fastq_file="$1"
output_fasta_file="${input_fastq_file%.fq}.fa"
awk '{if (NR%4==1) {sub(/^@/, ">", $0); print} else if (NR%4==2) print}' "$input_fastq_file" > "$output_fasta_file"
# 打印新 fasta 文件的开头
echo "The beginning of the new fasta file:"
head "$output_fasta_file"
# 检查行数，以免出错
echo "Check of line counts:"
fq_line_count=$(wc -l < "$input_fastq_file")
fa_line_count=$(awk 'NR%4==1 || NR%4==2' "$input_fastq_file" | wc -l)
half_fq_line_count=$((fq_line_count / 2))
echo "The original fastq file has $fq_line_count lines."
echo "The actual line count of the converted fasta file: $fa_line_count lines."
# 检查转换的行数是否约为原来的一半
if [ "$fa_line_count" -eq "$half_fq_line_count" ]; then
echo "Line count check passed."
else
echo "Line count check failed."
fi
运行示例：
➜ ./fq2fa SRR3091420_1_chr6.fastq
The beginning of the new fasta file:
>25
CTCAAGCAATCCTCCCACCTCAGCCTCCTGAGTAGTTGGGACTACAGGT
>1648
GGGATTGCAGGCATGCGCCACCATGCCCAGCTAATTTTTGTGTTTTTAG
>3092
CTGGAGTGCAATGGCGCGATCTCAGCTCACTGCAAGCTCTGCCTCCTGG
>3094
CACGCCACCACACCCAACTAATTTTTGTATTTTTAGTAGAGATGGGGTT
>9980
CCCCCGGTGTGTGATGTTCCCCTCCCTGTGTCCATGTGTTCTCATTGTT
Check of line counts:
The original fastq file has 3340676 lines.
The actual line count of the converted fasta file: 1670338 lines.
Line count check passed.
3. revCmpDNAseq
对DNA序列进行反转互补，用户可以选择仅反转或仅互补，或反转且互补；用户提供的可以是一条序列，或是一个文件（fa或fq），输出可以是fa/fq文件，请自行灵活地设计函数。
参数说明：
-r: 反转。对输入的DNA序列进行反转（例如：”GTCA” 变成 “ACTG”）。-c: 互补。对输入的DNA序列进行碱基互补（例如：”GTCA” 变成 “CAGT”）。-rc: 反转互补。先对序列进行互补，然后再反转（例如：”GTCA” 变成 “TGAC”）。-s <sequence>: 指定单条DNA序列。即可以直接在命令行中输入一个DNA序列进行处理。-f <file>: 指定输入文件（fasta 或 fastq 格式）。-o <output_file>: 指定输出文件。如果没有写-o或者-o后面没有指定输出文件名，还是会输出文件，输出的文件名默认为file_processed.filetype。
#!/bin/bash
# 打印帮助信息
usage() {
echo "Usage: $0 [-r] [-c] [-rc] [-s <sequence> | -f <file>] [-o <output_file>]"
echo "-r: Reverse the sequence."
echo "-c: Complement the sequence."
echo "-rc: Reverse complement the sequence."
echo "-s: Specify a single DNA sequence."
echo "-f: Specify an input file (fasta or fastq format)."
echo "-o: Specify an output file (optional)."
exit 1
}
# 反转序列
reverse_seq() {
echo "$1" | rev
}
# 互补序列
complement_seq() {
echo "$1" | tr 'ATCGatcg' 'TAGCtagc'
}
# 反转互补序列
reverse_complement_seq() {
echo "$1" | rev | tr 'ATCGatcg' 'TAGCtagc'
}
# 处理 fasta 或 fastq 文件
process_file() {
local file_path=$1
local operation=$2
local out_file=$3
while read -r line; do
if [[ $line == ">"* ]] || [[ $line == "@"* ]] || [[ $line == "+"* ]]; then
echo "$line"
else
$operation "$line"
fi
done < "$file_path" > "$out_file"
}
# 主程序
operation=""
input_seq=""
input_file=""
output_file=""
# 解析命令行参数
while getopts "rcs:f:o:" opt; do
case $opt in
r) operation="reverse_seq" ;;
c) operation="complement_seq" ;;
rc) operation="reverse_complement_seq" ;;
s) input_seq="$OPTARG" ;;
f) input_file="$OPTARG" ;;
o) output_file="$OPTARG" ;;
*) usage ;;
esac
done
# 校验参数
if [ -z "$operation" ] || { [ -z "$input_seq" ] && [ -z "$input_file" ]; }; then
usage
fi
# 如果提供了序列，直接处理序列
if [ -n "$input_seq" ]; then
if [ -n "$output_file" ]; then
$operation "$input_seq" > "$output_file"
else
$operation "$input_seq"
fi
# 否则处理文件
elif [ -n "$input_file" ]; then
if [ -z "$output_file" ]; then
output_file="${input_file%.*}_processed.${input_file##*.}"
fi
process_file "$input_file" "$operation" "$output_file"
fi
使用说明：
# 赋予可执行权限
chmod +x revCmpDNAseq
# 对单条序列进行操作
./revCmpDNAseq -r -s "GTCA"
# 对文件中的序列进行操作：
./revCmpDNAseq -rc -f input.fa -o output.fa
运行示例：
# 反转一条序列
➜ ./revCmpDNAseq -r -s "GTCAAAAAAATTTTTCGAA"
AAGCTTTTTAAAAAAACTG
# 反转互补一条序列
➜ ./revCmpDNAseq -rc -s "GTCAAAAAAATTTTTCGAA"
CAGTTTTTTTAAAAAGCTT
4. DNA2RNA
将fa或fq文件中的DNA序列转化成RNA序列。
#!/bin/bash
# 打印帮助信息
usage() {
echo "Usage: $0 [-s <sequence> | -f <file>] [output_file]"
echo "-s: Specify a single DNA sequence."
echo "-f: Specify an input file (fasta or fastq format)."
echo "Converts DNA sequences to RNA sequences."
exit 1
}
# 将DNA序列转换为RNA序列
convert_dna_to_rna() {
echo "$1" | tr 'Tt' 'Uu'
}
# 处理 fasta 或 fastq 文件
process_file() {
local file_path=$1
local out_file=$2
awk '{
if ($0 ~ /^>/ || $0 ~ /^@/ || $0 ~ /^\+/) {
print
} else {
gsub(/T/, "U");
gsub(/t/, "u");
print
}
}' "$file_path" > "$out_file"
}
# 主程序
input_seq=""
input_file=""
output_file="/dev/stdout"
# 解析命令行参数
while getopts "s:f:" opt; do
case $opt in
s) input_seq="$OPTARG" ;;
f) input_file="$OPTARG" ;;
*) usage ;;
esac
done
shift $((OPTIND -1))
# 检查额外的参数是否为输出文件
if [ "$#" -eq 1 ]; then
output_file="$1"
fi
# 校验参数
if { [ -n "$input_seq" ] && [ -n "$input_file" ]; }; then
usage
elif [ -z "$input_seq" ] && [ -z "$input_file" ]; then
usage
fi
# 处理输入
if [ -n "$input_seq" ]; then
convert_dna_to_rna "$input_seq" > "$output_file"
elif [ -n "$input_file" ]; then
process_file "$input_file" "$output_file"
fi
echo "Conversion complete. Output in $output_file"
使用说明：
# 赋予可执行权限
chmod +x DNA2RNA
# 处理单条序列：
./DNA2RNA -s "GTCA" output.fa
# 处理文件：
./DNA2RNA -f input.fa output.rna
脚本整合安装程序
参考gibbslab/biobash: Common BASH routines for computational biology (github.com)中的安装程序
#!/bin/bash
default="/usr/local/bin"
echo "
#
Bio_funclib INSTALLATION
#
"
failure=" == WARNING! Unable to install the tools =="
echo "Hit ENTER to use the default installation directory ($default) or provide a valid installation path below:"
read -p " > " dir
if [ -z "$dir" ]; then
installDir=$default
else
installDir=$dir
fi
# 检查目录是否存在
if [ -d "$installDir" ]; then
# 检查目录是否可写
if [ -w "$installDir" ]; then
echo "Installing to: $installDir"
# 复制文件到指定目录
cp ./DNA2RNA $installDir
cp ./fq2fa $installDir
cp ./revCmpDNAseq $installDir
cp ./statBioFile $installDir
# 确保所有脚本的执行权限
chmod +x ${installDir}/DNA2RNA
chmod +x ${installDir}/fq2fa
chmod +x ${installDir}/revCmpDNAseq
chmod +x ${installDir}/statBioFile
echo "Installation complete. Thank you for installing the tools."
else
echo $failure
echo "$installDir is not writable. Please make sure you have write access."
fi
else
echo $failure
echo "The directory you just typed ($installDir) does not exist."
echo "Please make sure route is correct"
fi
使用说明：
# 给予安装脚本执行权限
chmod +x Bio_funclib_install.sh
# 运行安装脚本
./install_script.sh
# 使用示例
DNA2RNA --help