格式转换

#对文件夹下的所有sra文件批量处理
for i in *sra
do
echo $i
# 对于双端测序数据，需要加--split-3参数，每样本处理约10min
fastq-dump --split-3 $i --gzip 
done
##fastaq转fasta
sed -n '1~4s/^@/>/p;2~4p' in.fastq > out.fasta

fastq格式解读

• 第一行以@开头，之后为序列的标识符以及描述信息（与FASTA格式的描述行类似）
• 第二行为序列信息
• 第三行以+开头，之后可以再次加上序列的标识及描述信息（可选）
• 第四行为质量得分信息，与第二行的序列相对应，长度必须与第二行相同

二代数据质检

for id in *fastq
do
echo $id
# 此处使用8线程，平均每文件处理约10min
fastqc -t 8 $id -o ./out
done

汇总质控结果

conda activate py37
multiqc -d qcout/ -o all/  #-d 是fastqc的结果文件夹 -o是结果文件夹

去除接头

for id in *fastq
do
echo $id
# 切除序列5'前10个碱基，每个文件处理约5min
cutadapt -u 10 -o ./cut_$id $id
done

表达量

for i in ./*.bam; do
    sorted_bam="${i%.bam}_sorted.bam"

    # Sort the BAM file
   samtools sort $i -o $sorted_bam

    # Index the sorted BAM file
    samtools index $sorted_bam
htseq-count -r pos -i Parent -t mRNA -m union -a 10 -f bam $sorted_bam ../Ppatens_318_v3.3.gene.gff3 >  $i.count

count_dir="./"
output_file="merged_htseq_counts.txt"

# 获取所有 .count 文件的列表
count_files=(${count_dir}/*.count)

# 使用第一个文件的基因名作为第一列，并开始构建结果文件
cut -f1 "${count_files[0]}" > gene_names.txt

# 生成合并文件的初始版本，包括基因名列
paste gene_names.txt <(cut -f2 "${count_files[0]}") > temp_merged.txt

# 遍历剩下的文件，将它们的计数值依次合并
for file in "${count_files[@]:1}"; do
  paste temp_merged.txt <(cut -f2 "$file") > temp_merged_new.txt
  mv temp_merged_new.txt temp_merged.txt
done

# 获取文件名并去掉 .count 后缀，作为列名
echo -e "Gene\\t$(basename -s .count "${count_files[@]}" | tr '\\n' '\\t' | sed 's/\\t$//')" > $output_file

# 将合并后的计数数据添加到输出文件中
cat temp_merged.txt >> $output_file

# 清理临时文件
rm gene_names.txt temp_merged.txt