Linux 生信学习路线¶

# === 目录操作 ===
pwd                          # 显示当前目录（你在哪里？）
ls                           # 列出当前目录的文件
ls -la                       # 详细列出所有文件（包括隐藏文件，显示大小和权限）
ls -lh                       # 文件大小用人能读懂的格式（KB/MB/GB）
cd /home/user/data           # 切换到指定目录
cd ..                        # 返回上一级目录
cd ~                         # 回到家目录
cd -                         # 回到上一次去过的目录

# === 创建和删除 ===
mkdir -p project/{data,results,scripts}  # 一次创建嵌套目录结构
# -p 表示自动创建父目录，{} 是bash展开语法
touch notes.txt              # 创建空文件（或更新文件时间戳）
cp file.txt backup/          # 复制文件
cp -r dir1/ dir2/            # 复制整个目录（-r = recursive 递归）
mv old_name.txt new_name.txt # 重命名文件
mv file.txt target_dir/     # 移动文件
rm file.txt                  # 删除文件（没有回收站！删了就没了！）
rm -r directory/             # 删除目录及其内容（⚠️ 谨慎使用）

# === 查看文件内容 ===
cat gene_list.txt            # 一次性输出全部内容（文件小时用）
head -20 expression.tsv      # 查看前20行（默认10行）
tail -20 expression.tsv      # 查看后20行
less huge_file.fastq         # 分页查看（按q退出，按/搜索）
wc -l reads.fastq            # 数行数（FASTQ行数÷4 = reads数）

# === 文本搜索（grep，生信天天用！）===
grep "BRCA1" gene_annotation.txt          # 搜索包含BRCA1的行
grep -c "BRCA1" gene_annotation.txt       # 只显示匹配的行数（-c = count）
grep -i "brca1" gene_annotation.txt       # 忽略大小写搜索（-i = ignore case）
grep -v "^#" variants.vcf                 # 排除以#开头的行（-v = invert，取反）
grep -w "TP53" gene_list.txt              # 精确匹配整个单词（-w = word）

# === 文件搜索 ===
find . -name "*.fastq.gz"                 # 在当前目录递归搜索所有fastq.gz文件
find . -name "*.bam" -size +1G            # 搜索大于1GB的BAM文件
find . -name "*.log" -mtime +7            # 搜索7天前修改的log文件

# === cut —— 按列提取（白话：从表格里切几列出来）===
cut -f1,3 gene_annotation.tsv             # 提取第1列和第3列（-f = field）
cut -f1-5 expression.tsv                  # 提取第1到5列
cut -d',' -f2 sample_info.csv             # 用逗号分隔（-d = delimiter），取第2列

# === sort —— 排序 ===
sort -k2,2nr expression.tsv              # 按第2列数值从大到小排
# -k2,2 表示按第2列排序
# n 表示按数值排（不是字符串）
# r 表示逆序（从大到小）

sort -k1,1 -k2,2n genes.bed              # 先按第1列字母排，再按第2列数值排

# === uniq —— 去重（必须先sort！）===
sort gene_list.txt | uniq                 # 去重
sort gene_list.txt | uniq -c              # 去重并计数（每个出现了几次）
sort gene_list.txt | uniq -d              # 只显示重复的行

# === awk —— 文本处理神器（面试必考！）===
# 基本语法：awk '{条件} {动作}' 文件

# 打印第1列和第3列
awk '{print $1, $3}' gene_annotation.tsv

# 筛选第5列（p值）小于0.05的行
awk '$5 < 0.05' results.tsv

# 计算第2列的平均值
awk '{sum += $2; n++} END {print sum/n}' expression.tsv

# 统计FASTQ文件的reads数（每4行一个read）
awk 'END {print NR/4}' reads.fastq

# 提取BED文件中chr17的区域，并计算区域长度
awk '$1 == "chr17" {print $0, $3-$2}' regions.bed

# === sed —— 流编辑器（替换文本）===
sed 's/old/new/g' file.txt               # 把所有old替换成new（g = global全局）
sed -n '10,20p' file.txt                 # 只打印第10到20行（-n + p = 只打印匹配的）
sed '/^#/d' file.vcf                     # 删除以#开头的行（d = delete）
sed 's/\t/,/g' data.tsv > data.csv       # tab转逗号（TSV转CSV）

# === 管道 | —— 把上一个命令的输出作为下一个命令的输入 ===

# 例1：统计VCF文件中有多少个变异（排除注释行）
grep -v "^#" variants.vcf | wc -l

# 例2：找出表达量最高的前10个基因
sort -k2,2nr expression.tsv | head -10

# 例3：统计每条染色体上有多少基因
cut -f1 gene_annotation.tsv | sort | uniq -c | sort -nr

# 例4：从FASTQ提取序列名并去重
grep "^@" reads.fastq | cut -d' ' -f1 | sort -u | wc -l

# === 重定向 ===
command > output.txt          # 输出写入文件（覆盖）
command >> output.txt         # 输出追加到文件
command 2> error.log          # 错误信息写入文件
command > output.txt 2>&1     # 标准输出和错误都写入同一个文件
command < input.txt           # 从文件读取输入

#!/bin/bash
# 脚本名：batch_fastqc.sh
# 功能：对所有FASTQ文件运行FastQC质控
# 用法：bash batch_fastqc.sh /path/to/fastq_dir /path/to/output_dir

set -euo pipefail    # 三个安全选项：
# -e: 任何命令出错就立即退出
# -u: 使用未定义变量时报错
# -o pipefail: 管道中任何命令出错就报错

# === 接收命令行参数 ===
INPUT_DIR=${1:?"错误：请提供FASTQ文件目录"}   # 第1个参数，没提供就报错
OUTPUT_DIR=${2:?"错误：请提供输出目录"}         # 第2个参数

# === 创建输出目录 ===
mkdir -p "${OUTPUT_DIR}"

# === 处理所有FASTQ文件 ===
for fq in "${INPUT_DIR}"/*.fastq.gz; do        # 遍历所有fastq.gz文件
    sample=$(basename "${fq}" .fastq.gz)        # 提取样本名（去掉路径和后缀）
    echo "[$(date '+%Y-%m-%d %H:%M:%S')] 正在处理: ${sample}"  # 打印时间和进度
    fastqc "${fq}" -o "${OUTPUT_DIR}" -t 4      # 运行FastQC，4线程
done

echo "全部完成！结果在 ${OUTPUT_DIR}"

#!/bin/bash
# 脚本名：rnaseq_pipeline.sh
# 功能：简化版RNA-seq分析流程（质控->比对->计数）

set -euo pipefail

# === 配置参数 ===
THREADS=8                                    # CPU线程数
GENOME_INDEX="/ref/STAR_index"               # STAR基因组索引
GTF="/ref/genes.gtf"                         # 基因注释文件
SAMPLES=("WT_1" "WT_2" "KO_1" "KO_2")       # 样本列表

# === 日志函数 ===
log() {
    echo "[$(date '+%H:%M:%S')] $*"          # 打印带时间戳的消息
}

# === 主流程 ===
for sample in "${SAMPLES[@]}"; do            # 遍历每个样本
    log "===== 开始处理样本: ${sample} ====="

    # 步骤1: 质控
    log "步骤1: FastQC 质控"
    fastqc "data/${sample}_R1.fastq.gz" "data/${sample}_R2.fastq.gz" \
        -o results/fastqc/ -t "${THREADS}"

    # 步骤2: 去接头（Trimming）
    log "步骤2: fastp 去接头和低质量"
    fastp \
        -i "data/${sample}_R1.fastq.gz" \    # 输入R1
        -I "data/${sample}_R2.fastq.gz" \    # 输入R2
        -o "results/trimmed/${sample}_R1.fq.gz" \  # 输出R1
        -O "results/trimmed/${sample}_R2.fq.gz" \  # 输出R2
        -h "results/fastp/${sample}.html" \  # HTML报告
        -w "${THREADS}"                       # 线程数

    # 步骤3: STAR比对
    log "步骤3: STAR 比对到参考基因组"
    STAR \
        --runThreadN "${THREADS}" \            # 线程数
        --genomeDir "${GENOME_INDEX}" \        # 基因组索引目录
        --readFilesIn "results/trimmed/${sample}_R1.fq.gz" \
                      "results/trimmed/${sample}_R2.fq.gz" \
        --readFilesCommand zcat \              # 输入是gz压缩的
        --outSAMtype BAM SortedByCoordinate \  # 输出排序后的BAM
        --outFileNamePrefix "results/star/${sample}_" \
        --quantMode GeneCounts                 # 同时输出基因计数

    # 步骤4: samtools 索引BAM
    log "步骤4: 索引BAM文件"
    samtools index "results/star/${sample}_Aligned.sortedByCoord.out.bam"

    log "===== 样本 ${sample} 处理完成 ====="
done

# 步骤5: MultiQC 汇总所有质控报告
log "步骤5: MultiQC 汇总报告"
multiqc results/ -o results/multiqc/

log "全部流程完成！"

# === 变量和字符串操作 ===
filename="sample_WT_1_R1.fastq.gz"
echo "${filename%.fastq.gz}"          # 去掉后缀：sample_WT_1_R1
echo "${filename##*/}"                # 去掉路径（取文件名）
echo "${filename/WT/KO}"              # 替换：sample_KO_1_R1.fastq.gz

# === 条件判断 ===
if [[ -f "results.tsv" ]]; then       # -f 检查文件是否存在
    echo "结果文件已存在"
elif [[ -d "results/" ]]; then        # -d 检查目录是否存在
    echo "结果目录存在但无文件"
else
    echo "需要先运行分析"
fi

# === 循环的几种写法 ===
# 写法1：for循环处理文件列表
for f in data/*.fastq.gz; do
    echo "处理: $f"
done

# 写法2：while循环读取文件每一行
while IFS=$'\t' read -r gene_id fc pval; do   # 按tab读取3列
    if (( $(echo "$pval < 0.05" | bc -l) )); then
        echo "${gene_id} 是显著差异基因"
    fi
done < results.tsv

# 写法3：并行处理（用 & 后台运行 + wait 等待）
for sample in WT_1 WT_2 KO_1 KO_2; do
    fastqc "data/${sample}.fastq.gz" -o results/fastqc/ &  # & 放后台
done
wait   # 等所有后台任务完成
echo "全部FastQC完成"

# === 安装 Miniforge（推荐，比 Anaconda 更适合生信）===
# wget https://github.com/conda-forge/miniforge/releases/latest/download/Miniforge3-Linux-x86_64.sh
# bash Miniforge3-Linux-x86_64.sh

# === 配置频道（channels）===
conda config --add channels defaults         # 默认频道
conda config --add channels bioconda          # 生信工具频道（必加！）
conda config --add channels conda-forge       # 社区频道（必加！）
conda config --set channel_priority strict    # 严格频道优先级

# === 创建环境 ===
# 推荐：每个项目/流程一个独立环境
conda create -n rnaseq python=3.10           # 创建名为rnaseq的环境，指定Python版本
conda activate rnaseq                         # 激活环境
conda deactivate                              # 退出环境

# === 安装工具 ===
# 用 mamba 代替 conda 安装（速度快10倍！）
mamba install -c bioconda fastqc multiqc star salmon   # 安装RNA-seq工具
mamba install -c bioconda samtools bedtools             # 安装常用工具

# === 环境管理 ===
conda env list                               # 列出所有环境
conda list                                   # 列出当前环境的所有包
conda env export > environment.yml           # 导出环境配置（可重复部署！）
conda env create -f environment.yml          # 从配置文件创建环境

# === 常见问题 ===
# 问题：安装冲突（dependency conflict）
# 解决：用mamba代替conda，或创建新的独立环境
mamba create -n variant_calling -c bioconda gatk4 picard samtools

# === 推荐的环境划分 ===
# 环境1：RNA-seq 分析
mamba create -n rnaseq -c bioconda \
    fastqc multiqc fastp star salmon featurecounts

# 环境2：宏基因组分析
mamba create -n metagenome -c bioconda \
    kneaddata metaphlan humann kraken2 bracken

# 环境3：机器学习
mamba create -n ml python=3.10 \
    scikit-learn pandas numpy matplotlib seaborn xgboost

# 环境4：R 分析
mamba create -n r_env -c conda-forge r-base=4.3 \
    r-tidyverse bioconductor-deseq2 bioconductor-edger

# === 记住：不同环境的工具不冲突 ===
conda activate rnaseq       # 做RNA-seq时激活这个
conda activate metagenome   # 做宏基因组时激活这个

#!/bin/bash
#SBATCH --job-name=star_align        # 任务名称
#SBATCH --partition=normal           # 队列/分区名
#SBATCH --nodes=1                    # 使用1个节点
#SBATCH --ntasks=1                   # 1个任务
#SBATCH --cpus-per-task=8            # 每个任务用8个CPU核心
#SBATCH --mem=32G                    # 申请32GB内存
#SBATCH --time=04:00:00              # 最大运行时间4小时
#SBATCH --output=logs/%x_%j.out      # 标准输出文件（%x=任务名,%j=任务ID）
#SBATCH --error=logs/%x_%j.err       # 错误输出文件

# === 加载环境 ===
source activate rnaseq               # 激活conda环境

# === 运行分析 ===
STAR \
    --runThreadN ${SLURM_CPUS_PER_TASK} \   # 使用申请的CPU数
    --genomeDir /ref/STAR_index \
    --readFilesIn data/sample_R1.fq.gz data/sample_R2.fq.gz \
    --readFilesCommand zcat \
    --outSAMtype BAM SortedByCoordinate \
    --outFileNamePrefix results/star/sample_

# === SLURM 常用命令 ===
sbatch job_script.sh                 # 提交任务
squeue -u $USER                      # 查看自己的任务队列
scancel 12345                        # 取消任务（12345是任务ID）
scancel -u $USER                     # 取消自己的所有任务
sinfo                                # 查看集群状态
sacct -j 12345 --format=JobID,State,Elapsed,MaxRSS  # 查看任务资源使用

# === 批量提交（数组任务）===
# 一次提交多个样本的分析，每个样本独立运行
#!/bin/bash
#SBATCH --array=0-3                  # 数组索引0到3（4个样本）
#SBATCH --cpus-per-task=8
#SBATCH --mem=32G

SAMPLES=("WT_1" "WT_2" "KO_1" "KO_2")
SAMPLE=${SAMPLES[$SLURM_ARRAY_TASK_ID]}  # 用数组索引取样本名

echo "正在处理样本: ${SAMPLE}"
# ... 后续分析命令用 ${SAMPLE} 替代样本名

# Snakefile —— 用Python语法定义分析流程

# 配置
SAMPLES = ["WT_1", "WT_2", "KO_1", "KO_2"]

# 最终目标
rule all:
    input:
        "results/multiqc/multiqc_report.html",          # 质控报告
        expand("results/counts/{sample}.count", sample=SAMPLES)  # 所有计数文件

# 步骤1：FastQC质控
rule fastqc:
    input:
        r1 = "data/{sample}_R1.fastq.gz",               # 输入文件（用通配符）
        r2 = "data/{sample}_R2.fastq.gz"
    output:
        html1 = "results/fastqc/{sample}_R1_fastqc.html",
        html2 = "results/fastqc/{sample}_R2_fastqc.html"
    threads: 2                                           # 使用2个线程
    shell:
        "fastqc {input.r1} {input.r2} -o results/fastqc/ -t {threads}"

# 步骤2：STAR比对
rule star_align:
    input:
        r1 = "data/{sample}_R1.fastq.gz",
        r2 = "data/{sample}_R2.fastq.gz"
    output:
        bam = "results/star/{sample}_Aligned.sortedByCoord.out.bam"
    params:
        index = "/ref/STAR_index",                       # 参数（不算输入文件）
        prefix = "results/star/{sample}_"
    threads: 8
    shell:
        """
        STAR --runThreadN {threads} \
             --genomeDir {params.index} \
             --readFilesIn {input.r1} {input.r2} \
             --readFilesCommand zcat \
             --outSAMtype BAM SortedByCoordinate \
             --outFileNamePrefix {params.prefix}
        """

# 步骤3：featureCounts计数
rule featurecounts:
    input:
        bam = "results/star/{sample}_Aligned.sortedByCoord.out.bam"
    output:
        count = "results/counts/{sample}.count"
    params:
        gtf = "/ref/genes.gtf"
    threads: 4
    shell:
        "featureCounts -T {threads} -a {params.gtf} -o {output.count} {input.bam}"

# === Snakemake 运行命令 ===
snakemake -n                         # 干跑（dry-run），查看会执行什么
snakemake -j 8                       # 运行，最多同时用8个核心
snakemake --dag | dot -Tpdf > dag.pdf  # 生成流程图
snakemake --use-conda                # 自动管理conda环境
snakemake --cluster "sbatch -c {threads} --mem={resources.mem}"  # 提交到集群