标签：分析 samtools name -- chip file bam id

参考博客

https://zhuanlan.zhihu.com/p/377600056
https://www.jianshu.com/p/96688fecd864
https://zhuanlan.zhihu.com/p/676395563

查看质控

# 使用fastqc查看质控结果
fastqc -t 40 -o 1_rawdata_qc/ 0_rawdata/*.fastq.gz 1> 1_fastqc.log

# 使用MutliQC整合FastQC结果
# multiqc 1_rawdata_qc/

使用trim_galore软件过滤低质量的reads

ls 0_rawdata/*.fastq.gz | while read id;
do 
nohup trim_galore -q 25 --phred33 --length 25 -e 0.1 --stringency 4 -o 2_cleandata $id & 
done

查看过滤后的质控

# 使用fastqc查看质控结果
fastqc -t 40 -o 3_cleaneddata_qc/ 2_cleandata/*.fq.gz 1> 3_fastqc.log

# 使用MutliQC整合FastQC结果
# multiqc 3_cleaneddata_qc/

使用Bowtie2软件比对reads到参考基因组

# 软链接参考基因组和索引文件
ln -s /vol1/yudawei/cuttag/genome/pig/pig.fa* /vol1/yudawei/chip/genome

# 比对
ls 2_cleandata/*_1_trimmed.fq.gz | while read id;
do
file=$(basename $id) # 使用 basename 命令提取变量 id 中的文件名（去除路径部分）
sample=${file%%_trimmed.*} # 从变量 file 中提取样本名，移除文件名中 _trimmed 及其后面的部分
bowtie2 -p 4 -x /vol1/yudawei/chip/genome/pig -1 2_cleandata/$sample_1_trimmed.fq.gz -2 2_cleandata/$sample_2_trimmed.fq.gz | samtools sort -O bam -@ 20 -o - > 3_bam/${sample}.bam
done

对bam文件进行质控

ls 3_bam/*.bam | xargs -i samtools index {} 
ls 3_bam/*.bam | while read id ;
do 
nohup samtools flagstat $id > $(basename $id ".bam").stat & 
done

去除PCR重复

先将sam文件转化为bam文件，然后进行排序(按照坐标和nam)，利用排好序的bam文件来进行fixmate校正以及去重，双端测序数据用samtools rmdup效果很差，相对来说更建议用picard工具的MarkDuplicates功能。samtools的rmdup是直接将这些重复序列从比对BAM文件中删除掉，而Picard的MarkDuplicates默认情况则只是在BAM的FLAG信息中标记出来，而不是删除，因此这些重复序列依然会被留在文件中，只是我们可以在变异检测的时候识别到它们，并进行忽略。虽然samtools过时了，但不过好在可以使用samtools markdup代替。samtools markdup与picard MarkDuplicates采用类似的去重策略。

dir="bowtie2_output"
cd "$dir"
for sam_file in *.sam
do
    echo "$(date): Starting processing of $sam_file"

    name="${sam_file%.sam}"

    samtools view -bS "$sam_file" -o "$name.bam" -@ 8
    samtools sort -n "$name.bam" -o "$name.namesrt.bam" -@ 8
    samtools fixmate -m "$name.namesrt.bam" "$name.fixmate.bam" -@ 8
    samtools sort "$name.fixmate.bam" -o "$name.coordsrt.bam" -@ 8
    samtools markdup -r "$name.coordsrt.bam" "$name.rmdup.bam" -@ 8

    rm "$name.bam" "$name.namesrt.bam" "$name.fixmate.bam"
    echo "$(date): Finished processing of $sam_file"
done
echo "All files have been processed."

观察去重和未去重的差异

dir="bowtie2_output"
cd "$dir"

for sorted_file in *coordsrt.bam
do
    base="${sorted_file%.coordsrt.bam}"
    rmdup_file="${base}.rmdup.bam"
    samtools flagstat "${sorted_file}" > "${sorted_file}.flagstat"
    samtools flagstat "${rmdup_file}" > "${rmdup_file}.flagstat"

    echo "Comparing flagstat reports for ${base}" >> comparison.log
    diff -y "${sorted_file}.flagstat" "${rmdup_file}.flagstat" >> comparison.log
done
echo "Flagstat comparison is complete."

bamCompare

为了能够比较不同的样本，需要对先将基因组分成等宽分箱(bin)，统计每个分箱的read数，最后得到描述性统计值

bamCompare -p 10 -b1 <sample>.rmdup.bam -b2 <sample>._IN.rmdup.bam --skipNAs --scaleFactorsMethod readCount --operation subtract --outFileFormat bedgraph -o ../bw_data/<sample>.compareInput.bedgraph --extendReads 200
# -p为线程数
# -b1为ip
# -b2为input
# --skipNA是跳过bin中没有覆盖的部分
# --scaleFactorsMethod是选择缩放样本的方法

# 可能会存在负值，所以这一步需要将其矫正为0
awk '{if($4<0){$4=0};print}' <sample>.compareInput.bedgraph > <sample>.move0.bedgraph

totalreads=`awk '{sum+=$4}END{print sum}' <sample>.move0.bedgraph` |echo 'reads left after remove duplication: '$totalreads

awk -v totalreads=$totalreads '{$4=$4/totalreads*1000000;print}' <sample>.move0.bedgraph >  <sample>.rpm.bedgraph

sort -k1,1 -k2,2n  <sample>.rpm.bedgraph  >  <sample>.sort.bedgraph

使用macs2软件进行ChIP-seq peak calling

ls 4_dedup/*.markdup.bam | while read id
do 
nohup macs2 callpeak -f BAM -B -g 2.2e9 -n $id -p 0.05 --outdir 5/ >./$id.log &
done

标签：分析,samtools,name,--,chip,file,bam,id
From： https://www.cnblogs.com/cauwj/p/18137153

chip分析