- 查看文件夹大小：
du -sh
脑图
1. 数据准备
- 参考基因组: genome.fa
- 将参考基因组的格式改为">chr1"，而不是单纯的数字。
sed 's/>/>chr/' genome.fa
- gtf注释文件: genes.gtf
- 提取"protein_coding"
grep "protein_coding" genes.gtf > genes.filter.gtf
- 提取最长转录本（待后面讲）
怎么看有没有可变剪接？
awk '$3=="gene"' Danio_rerio.GRCz11.113.filter.gtf|wc
awk '$3=="transcript"' Danio_rerio.GRCz11.113.filter.gtf|wc
转录本数目远远多于基因数目，说明这个物种注释到了转录本的水平。
可以从基因层面研究，也可以从转录本层面研究。如果从基因层面研究，则TF注释到基因上就行了，不用管TF是位于基因哪个转录本的上游，为什么不用去管？一个基因对应多个转录本，多个转录本的位置在gtf文件中可能并不准确,尤其是对于非模式物种来说。
提取最长转录本代表这个基因，每个基因提取1个。
2. 软件安装
参见：tutorial的book文档，这里不多赘述。
#使用容器中的bowtie2
singularity exec bowtie2_2.5.2.sif bowtie2
3. 操作
- 医研院服务器的容器加载
module load singularity/3.8.0
- 加载已经打包好在容器中的软件
singularity exec bowtie2_2.5.2.sif bowtie2
测试一下流程
note：研究的物种是拟南芥。
mdkir testPeakSnake
cd testPeakSnake/
ln -s ../01raw_data #后面出错，说明要拷贝过来，不能链接
cp -r ../02ref/ ./
cp ../PeakSnake/Snakefile ./
cp ../PeakSnake/config.yaml ./
snakemake -np #假装跑一下这个流程，看有哪些步骤.
#真正跑
cp ../PeakSnake/run.sh ./
cat run.sh
nohup snakemake --cores 46 --use-singularity --keep-going &
htop -u wangmengj #查看与运行情况
#出错
cp -r ../01raw_data ./ #解决，因为是容器，所有数据必须位于当前目录下。
1.第一步：测序数据过滤和质控
1.去除低质量reads: 如果1条reads有150bp，中间有10bp的低质量。
两种选择：①整条都丢掉即remove；②只丢掉中间10bp即trim。
大多数情况下，采用的是remove这种方式，因为现在测序价格较低。
具体用法见基因课tutorial。
Q：不指定adapter序列，软件怎样检测出来序列？不指定它咋知道adapter长什么样子？
先说单末端测序，原则上数据在下机的时候，加的接头序列就已经被去干净了。如果有一些没去干净，fastp软件会发现在序列末尾，有一些固定的碱基出现的频数很高，它就认为应该是adapter。
双末端测序就更简单了，双末端只有在DNA片段太短了，测通了，才会测到adapter。这种情况，只需要将read1和read2进行比较，重叠了，重叠的部分就是真正的序列。多出来的就是adapter。
note: 使用快捷命令，只在当前会话终端中有效，临时命令
alias le='less -SN'
le quality_control/fastp_summary.tsv
csvlook quality_control/fastp_summary.tsv | le
得到这样的结果：
2.第二步：比对到参考基因组
二代数据比对到参考基因组，最早用blast，比对比较长的序列，是局部比对，分区段的。后来，华大基因开发的SOAP。后来李恒开发的bwa，取代blast, SOAP，成为最主流的比对软件，现在十多年过去了，基因组重测序还是用bwa。但都不能将RNA测序数据比对到参考基因组，开gap的比对，因为RNA有内含子，需要跳过一大段去比对。
后来开发针对RNA-seq的软件：tophat，至少用了五六年时间，后来tophat作者又开发hisat，后来又升级一个版本叫hisat2。
后来有STAR软件，有个ENCODE计划，启动项目前对所有软件测评，测评完说STAR是最好的转录组比对软件。
像我们用的核糖体RNA, small RNA也是用bowtie2，因为small RNA很短,中间没有内含子。
- bowtie2的结果（单端数据）
Q：什么是index？
答：比对前，对参考基因组做的一些准备工作，比如：把染色体长度从长到短排序，计算每条染色体多长，这种准备工作统称index，类似黑箱子。
- sam文件格式
是一个文本格式。
ChIP-Seq在构建文库过程中，对起始量要求特别高，需要几百万个细胞，而且效率特别低，经常导致文库质量不高，比如: 数据质量不高，duplication比例太高等各种问题。
- dedup
Q：什么是duplicates?
答：基本上是由于PCR产生的。测序是一个随机抽样的过程。
picard MarkDuplicates，检测比对结果中，两条序列头也一样，尾也一样，长得完全一样，比对到同一个地方，一个碱基不差，这种就把它认定成duplicates。
picard MarkDuplicates I=raw_data/test.bam O=clean_datatest.rmdup.bam M=matrix.txt REMOVE_DUPLICATES=False #这个参数表示只标记，不去除
- 使用deeptools
module load singularity/3.8.0
cd /home/wangmengj/01_genek/01_ChIP-seq/PeakSnake/sifs
singularity exec deeptools_20231220.sif
singularity exec deeptools_20231220.sif alignmentSieve -h
- sam flag是什么？
谷歌搜索sam explain，进入broad institute
在sam文件里，每一行比对都有一个flag值，这个flag值描述了这个reads是如何比对到参考基因组上的，以及这个reads的状态，存在于sam文件的第二列。
- 过滤unmapped
- 过滤次级比对，比如一个reads比对到两个地方，但是第二个地方不是最好的地方。
以前很多文章里，会去除多处比对的reads，其实多处比对的reads在bam里面是不标记的，是需要自己去统计的。像这样简单去除一下不是最优的比对就可以了。
- 去除黑名单
加上线粒体和叶绿体：
看线粒体多长samtools faidx
samtools faidx genome.fa
less genome.fa.fai #第二列就是染色体长度
Q: 非模式物种的黑名单怎么办？
1.blacklist
2.greenscreen
igv可视化
- 有效基因组大小：比如基因组上有很多N，这种的话reads就比对不上去，就不算大小。
1.deeptools提供的
2.自己算
# 非要绝对路径
/home/wangmengj/01_genek/01_ChIP-seq/PeakSnake/scripts/faCount -summary ref/genome.fa
# 然后将总的数量减去N的数量
然后就可以导入bigwig查看即可。
- 设置track高度:设成120
- 设置data range:设成120
